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Abstract 


This  report  details  the  development  of  an  A-channel  spatial  power  spectrum  estimation 
technique  called  the  SPatial  /ncoherent  Region  Estimator  (SPIRE).  It  was  developed  in 
support  of  research  aimed  at  characterizing  the  spatial  spreading  of  HF  signals  caused 
by  the  high  latitude  propagation  environment  with  the  ultimate  goal  of  improving  high 
latitude  HF  radio  direction  finding  performance.  Based  on  the  assumption  of  temporal 
and  spatial  incoherence,  SPIRE  uses  a  surprisingly  simple,  yet  effective,  approach 
based  on  maximum  likelihood  principles  to  model  the  spatial  power  spectrum.  The 
result  is  an  algorithm  which  provides  a  more  accurate  and  more  informative 
characterization  of  the  spatial  nature  of  incoming  signals  than  currently  popular 
conventional  and  modem  superresolution  algorithms.  This  characterization  includes 
bearing,  spatial  extent  and  power  distribution,  and  total  power  of  the  signal.  It  also 
includes  total  noise  power  and  modeling  accuracy.  The  performance  and  utility  of  the 
SPIRE  algorithm  is  demonstrated  using  both  simulation  and  off-air  data. 


Resume 


Ce  rapport  decrit  l’elaboration  d’une  technique  d’estimation  du  spectre  de  puissance 
spatial  utilisant  N  canaux  appelee  SPIRE  (SPatial  Incoherent  Region  Estimator,  c.-a-d. 
estimateur  de  region  spatiale  incoherente).  Cette  technique  a  ete  elaboree  pour  appuyer 
la  recherche  visant  a  caracteriser  l’etalement  des  signaux  HF  cause  par  le  milieu  de 
propagation  aux  hautes  latitudes  dans  le  but  premier  d’ameliorer  les  performances  de  la 
radiogoniometrie  HF  a  ces  latitudes.  A  partir  de  l’hypothese  d’incoherence  temporelle 
et  spatiale,  SPIRE  utilise  une  approche  etonnamment  simple,  mais  efficace,  basee  sur 
les  principes  du  maximum  de  vraisemblance,  pour  modeliser  le  spectre  de  puissance 
spatial.  Le  resultat  est  un  algorithme  qui  permet  une  caracterisation  plus  precise  et  plus 
informative  de  la  nature  spatiale  des  signaux  re?us  que  les  algorithmes  classiques  et 
modemes  de  superresolution  qui  sont  actuellement  repandus.  Cette  caracterisation 
porte  sur  le  relevement,  sur  la  distribution  de  puissance  et  l’etendue  spatiales,  ainsi  que 
sur  la  puissance  totale  du  signal.  Elle  porte  aussi  sur  la  puissance  de  bruit  totale  et  la 
precision  de  la  modelisation.  Les  performances  et  l’utilite  de  l’algorithme  SPIRE  sont 
demontrees  par  simulation  et  par  utilisation  de  donnees  captees  en  direct. 
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Executive  summary 


The  requirement  exists  to  improve  the  accuracy  of  strategic  high  latitude  HF  direction 
finding  (DF)  systems.  In  the  past,  poor  DF  accuracy  derived  from  Arctic  measurements 
has  led  to  low  confidence  in  high  latitude  sites  despite  the  strategic  relevance  of  these 
sites  for  transmitter  geolocation. 

Disturbances  in  the  ionosphere  which  scatter  the  signal  over  a  range  of  azimuth  and 
elevation  directions  (spatial  spreading)  are  a  major  problem.  This  leads  to  situations 
which  are  not  well  modeled  by  currently  popular  DF  algorithms  since  the  signal  is 
assumed  to  be  cleanly  reflected  from  the  ionosphere.  The  result  can  be  bearings  with 
large  biases  and  standard  deviations. 

Recently  a  new  approach,  the  Spread  Maximum  Likelihood  (SML)  algorithm,  was 
developed,  incorporating  a  spatial  spreading  model.  Despite  improved  performance 
compared  to  other  DF  algorithms,  performance  was  still  not  as  good  as  desired.  The 
problem  may  be  that  the  SML  algorithm  employs  a  signal  scattering  model  requiring 
the  spatial  shape  (but  not  size)  and  power  density  profile  of  the  scattering  region  to  be 
constant  and  known  a  priori.  Not  only  has  the  optimum  shape  not  yet  been  determined, 
there  is  reason  to  believe  that,  given  the  dynamic  nature  of  the  high  latitude  ionosphere, 
the  shape  of  the  scattering  region  may  change  significantly  over  time. 

To  understand  the  true  nature  of  these  scattering  regions  better,  the  SML  algorithm  has 
been  radically  modified  to  include  estimation  of  the  spatial  shape  and  power  density  of 
the  received  signal.  The  result  is  a  new  algorithm  called  the  SPatial  /ncoherent  Region 
Estimator  (SPIRE),  which  effectively  maps  out  the  spatial  power  spectrum  of  the  radio 
sky.  Surprisingly,  despite  the  increase  in  number  of  signal  parameters  to  be  estimated, 
the  SPIRE  algorithm  is  much  simpler  and  faster  than  the  SML  algorithm.  In  fact,  the 
SPIRE  algorithm  is  comparable  in  speed  to  the  current  generation  of  superresolution 
algorithms  and  is  much  more  suited  to  the  problem  of  spatially  spread  signals. 

The  development  and  evaluation  of  the  performance  of  the  SPIRE  algorithm  is 
documented  in  this  report.  The  performance  evaluation  was  carried  out  using  both 
simulated  and  off-air  data  and  also  included  comparisons  with  other  superresolution  DF 
algorithms. 

In  comparisons  using  simulated  data,  the  effects  of  noise,  spatial  spreading,  and  the 
shape  of  the  spread  region  were  considered,  as  well  as  the  ability  to  detect  a  weaker 
point-source  (no  spatial  spreading)  signal  in  the  presence  of  a  stronger  spatially  spread 
signal.  For  the  signal  conditions  tested,  the  performance  of  the  SPIRE  algorithm  was 
generally  found  to  be  as  good  as,  or  better  than,  the  performances  of  currently  popular 
superresolution  DF  algorithms. 

In  comparisons  using  off-air  data  collected  with  the  Vortex  system  (a  multi-channel  HF 
receiver  system  located  at  CFS  Alert)  in  1995-96,  the  results  were  inconclusive. 
However,  the  usefulness  of  the  SPIRE  algorithm  to  detect  problems,  particularly  the 
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ability  to  measure  modeling  error,  made  it  possible  to  show  that  antenna  mutual 
coupling  problems  had  almost  certainly  corrupted  the  measurements.  Ideally,  mutual 
coupling  effects  could  be  corrected  by  doing  advanced  processing  to  calibrate  the  data; 
however,  more  research  work  needs  to  be  done  in  this  area. 


W.J.L.Read.  2000.  A  Spatial  Power  Spectrum  Estimator  For  Distributed  Signals.  DREO  TR 
2000-099.  Defence  Research  Establishment  Ottawa. 
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Sommaire 


II  est  necessaire  d’ameliorer  la  precision  des  systemes  strategiques  de  radiogoniometrie 
HF  aux  hautes  latitudes.  Par  le  passe,  la  faible  precision  obtenue  en  radiogoniometrie 
avec  les  mesures  faites  dans  l’Arctique  a  mene  a  un  bas  niveau  de  confiance  dans  les 
sites  des  hautes  latitudes,  malgre  la  pertinence  strategique  de  ces  sites  pour  etablir  la 
position  geographique  des  emetteurs. 

Un  probleme  majeur  est  lie  aux  perturbations  ionospheriques  qui  diffusent  le  signal  sur 
une  gamme  de  directions  en  azimut  et  en  site  (etalement  spatial).  II  en  resulte  des 
situations  qui  ne  sont  pas  bien  modelisees  par  les  algorithmes  de  radiogoniometrie 
actuellement  repandus,  etant  donne  que  le  signal  est  considere  par  hypothese  comme 
etant  nettement  reflechi  par  1’ ionosphere.  Les  relevements  obtenus  peuvent  alors 
comporter  des  biais  et  des  ecarts-types  importants. 

Recemment,  un  nouvel  outil,  1’algorithme  a  vraisemblance  maximale  etalee  (SML), 
comprenant  un  modele  d’etalement  spatial,  a  ete  elabore.  Bien  que  cet  algorithme  soit 
plus  efficace  que  les  autres  algorithmes  de  radiogoniometrie,  l’efficacite  desiree  n’etait 
pas  encore  atteinte.  Le  probleme  peut  Stre  du  au  fait  que  1’ algorithme  SML  utilise  un 
modele  de  diffusion  du  signal  dans  lequel  la  forme  spatiale  (mais  pas  la  taille)  et  le 
profil  de  densite  de  puissance  de  la  region  de  diffusion  doivent  etre  constants  et  connus 
a  priori.  Non  seulement  la  forme  optimale  n’a  pas  encore  ete  determinee,  mais  on  est 
justifie  de  croire  qu’en  raison  de  la  nature  dynamique  de  l’ionosphere  aux  hautes 
latitudes,  la  forme  de  la  region  de  diffusion  peut  changer  beaucoup  avec  le  temps. 

Pour  mieux  comprendre  la  nature  reelle  de  ces  regions  de  diffusion,  on  a  modifie 
radicalement  1’algorithme  SML  de  fagon  a  y  inclure  une  estimation  de  la  forme  spatiale 
et  de  la  densite  de  puissance  du  signal  re<ju,  ce  qui  a  donne  un  nouvel  algorithme  appele 
SPIRE  (SPatial  Incoherent  Region  Estimator,  c.-a-d.  estimateur  de  region  spatiale 
incoherente),  qui  represente  efficacement  le  spectre  de  puissance  spatial  de  la 
propagation  radio  ionospherique.  Etonnamment,  en  depit  de  l’accroissement  du  nombre 
de  parametres  de  signaux  a  evaluer,  1’algorithme  resultant  est  beaucoup  plus  simple  et 
plus  rapide  que  1’ algorithme  SML,  avec  une  vitesse  comparable  a  celle  des  algorithmes 
de  superresolution  de  la  generation  actuelle,  mais  beaucoup  mieux  adapte  au  probleme 
des  signaux  a  etalement  spatial. 

L’ elaboration  de  1’algorithme  SPIRE  est  documentee  dans  ce  rapport  et  ses 
performances  sont  evaluees  a  l’aide  de  donnees  simulees.  Cette  evaluation  tient  compte 
des  effets  du  bruit,  du  degre  d’etalement  spatial  et  de  la  forme  de  la  region  d’etalement. 
Elle  permet  aussi  de  detecter  un  signal  de  source  ponctuelle  (pas  d’etalement  spatial) 
faible  en  presence  de  signaux  a  etalement  spatial  plus  forts. 

Dans  les  comparaisons  effectuees  a  l’aide  de  donnees  simulees,  les  performances  de 
1’algorithme  SPIRE  se  sont  generalement  revelees  egales  ou  superieures  a  celles  des 
algorithmes  de  radiogoniometrie  actuellement  repandus,  pour  les  diverses  conditions  de 
signaux  ayant  fait  l’objet  d’essais. 
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Dans  les  comparaisons  effectuees  a  l’aide  de  donnees  en  direct  recueillies  avec  le 
sy steme  de  reception  HF  multi voie  Vortex  a  la  SFC  Alert  en  1995-1996,  les  resultats 
n’etaient  pas  concluants.  Cependant,  l’utilite  de  l’algorithme  SPIRE  pour  detecter  les 
problemes,  en  particular  la  capacite  de  mesurer  l’erreur  de  moderation,  a  permis  de 
montrer  que  les  problemes  de  couplage  mutuel  d’antenne  avaient  presque  certainement 
corrompu  les  mesures.  Idealement,  les  effets  du  couplage  mutuel  pourraient  etre 
corriges  par  un  traitement  evolue  servant  a  etalonner  les  donnees,  mais  d’autres  travaux 
de  recherche  sont  necessaires  dans  ce  domaine. 


WJ.L.Read.  2000.  Une  technique  d’estimation  du  spectre  de  puissance  spatial  pour  les 
signaux  distribues.  DREO  TR  2000-099.  Centre  de  recherches  pour  la  defense,  Ottawa. 
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1.  INTRODUCTION 


The  requirement  exists  to  improve  the  accuracy  of  strategic  HF  direction  finding 
systems,  particularly  in  the  Arctic.  In  the  past,  poor  direction  finding  (DF)  accuracy 
derived  from  Arctic  measurements  has  led  to  low  confidence  in  high  latitude  sites 
despite  the  strategic  relevance  of  these  sites  for  transmitter  geolocation. 

Patches  of  enhanced  electron  density  and  associated  instabilities  in  the  F  layer  of  the 
ionosphere,  which  drift  across  the  polar  cap  during  darkness  in  a  roughly  antisunward 
direction  at  speeds  ranging  from  a  few  hundreds  to  over  one  thousand  meters  per 
second  [1],  are  a  major  problem.  These  patches  can  cause  scattering  of  a  signal  from 
azimuth  directions  which  are  very  different  from  the  true  bearing  of  the  transmitter. 

One  avenue  of  investigation  being  pursued  is  the  development  of  new  DF  algorithms 
which  are  better  matched  to  the  high  latitude  HF  signal  environment.  Currently  popular 
DF  algorithms  assume  that  the  incoming  signal  can  be  modeled  as  a  point  source,  or 
equivalently,  that  the  incoming  radio  signal  has  a  planar  wavefront.  This  is  reasonable 
if  the  size  of  the  transmission  source  is  extremely  small  relative  to  its  range,  the  size  of 
the  DF  array  is  also  small  relative  to  the  transmitter  range,  the  ionosphere  acts  as  a 
perfect  or  near-perfect  reflector,  and  local  site  multipath  can  be  ignored.  Unfortunately, 
at  high  latitudes,  during  periods  when  scattering  from  large  moving  patches  occurs,  the 
received  signal  arrives  from  a  range  of  bearings  in  both  azimuth  and  elevation.  In  this 
report,  signals  of  this  type  are  called  spread-source  signals.  Typical  modem  DF 
algorithms  estimate  a  spread-source  signal  as  being  a  cluster  of  several  point-source 
signals  coming  from  the  same  direction. 

Although  representation  of  a  spread-source  signal  as  a  cluster  of  several  point-source 
signals  is  useful,  the  information  provided  is  degraded.  For  example,  it  becomes  more 
difficult  to:  determine  the  spatial  extent  of  the  signal;  determine  if  two  or  more  similar 
bearing  estimates  represent  different  reflection  paths  or  a  single  scattering  region  of  the 
ionosphere;  and,  detect  and  estimate  the  direction  of  weaker  signals.  Solving  this  latter 
difficulty  is  important  since  a  previous  study  [2]  has  shown  that  sporadic-E  propagation, 
when  it  exists  and  can  be  measured,  leads  to  good  bearing  estimates.  Hence,  it  is 
important  that  a  DF  algorithm  be  able  to  detect  and  determine  the  bearing  of  a  weaker 
sporadic-E  reflected  signal  in  the  presence  of  one  or  more  stronger  F-region  signals. 

Recently  a  new  approach,  the  Spread  Maximum  Likelihood  (SML)  algorithm  [3],  was 
developed,  incorporating  a  spatial  spreading  model  to  handle  scattering  effects  due  to 
the  ionosphere.  Despite  improved  performance  compared  to  other  DF  algorithms  [1], 
[3],  [4],  performance  was  still  not  as  good  as  desired.  The  problem  may  be  that  this  new 
algorithm  employs  a  signal  scattering  model  requiring  the  shape  (but  not  size)  and 
power  density  profile  of  the  scattering  region  to  be  constant  and  known  a  priori.  Not 
only  has  the  optimum  shape  not  yet  been  determined,  there  is  reason  to  believe  that, 
given  the  dynamic  nature  of  the  high  latitude  ionosphere,  the  shape  of  the  scattering 
region  may  change  significantly  over  time. 
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To  understand  the  true  nature  of  these  scattering  regions  better,  the  SML  algorithm  has 
been  radically  modified  to  include  estimation  of  the  spatial  shape  and  power  density  of 
the  received  signal.  The  result  is  a  new  algorithm  called  the  SPatial  /ncoherent  Region 
Estimator  (SPIRE),  which  effectively  maps  out  the  spatial  power  spectrum  of  the  radio 
sky.  Surprisingly,  despite  the  increase  in  number  of  signal  parameters  to  be  estimated, 
the  SPIRE  algorithm  is  much  simpler  and  faster  than  the  SML  algorithm.  In  fact,  the 
SPIRE  algorithm  is  comparable  in  speed  to  the  current  generation  of  superresolution 
algorithms  and  is  much  more  suited  to  the  problem  of  spatially  spread  signals,  as  will 
be  shown  through  simulation. 

An  unexpected  bonus  of  the  SPIRE  algorithm  approach  is  a  function  which  effectively 
measures  the  modeling  error.  This  error  function  provides  a  useful  way  to  determine 
proper  sampling  size,  and  whether  local  multipath/coupling  effects  may  be  degrading 
the  results.  This  is  demonstrated  using  high  latitude  HF  data  collected  with  the  Vortex 
system  (an  experimental  12  channel  HF  collection  system)  in  1995-96  at  CFS  Alert. 

The  layout  of  the  rest  of  the  report  is  as  follows.  In  Section  2,  the  underlying  maximum 
likelihood  approach  is  discussed  and  the  appropriate  cost  and  model  error  functions  are 
introduced.  This  is  followed  by  the  development  of  a  signal  model  which  reduces  the 
estimation  problem  to  calculation  of  the  signal  and  noise  power  parameters  only.  In 
Section  3,  the  basic  algorithmic  procedure  for  signal  and  noise  power  parameter 
estimation  is  developed.  In  Section  4,  problems  with  spatial  ambiguities  are  dealt  with 
through  the  adoption  of  a  “simpler  is  better”  guideline.  An  implementation  scheme  is 
also  introduced  which  leads  to  significantly  faster  processing.  Additionally  issues,  such 
as  sample  size  and  maximum  number  of  signals,  are  also  addressed.  In  Section  5, 
simulated  data  is  used  to  compare  the  performance  of  the  SPIRE  algorithm  with  other 
popular  DF  algorithms.  In  Section  6,  performance  is  again  compared  using  off-air  data. 
This  data  is  also  processed  to  show  the  usefulness  of  measuring  modeling  error  for 
sample  size  determination  and  mutual  coupling/local  site  multipath  detection.  Finally, 
in  Section  7,  the  conclusions  and  recommendations  are  presented. 
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2.  DEVELOPMENT  OF  THE  SPIRE  ALGORITHM 

2.1  Maximum  Likelihood  Estimation 


A  successful,  albeit  often  computationally  intensive,  approach  to  estimation  is  based  on 
the  maximum  likelihood  method.  Essentially  the  idea  is  to  find  the  most  likely  state  of  a 
signal  process  given  a  set  of  measurement  observations  made  of  the  process.  Since  this 
is  a  statistical  approach,  the  method  applies  to  cases  where  the  signal  generation  is  a 
random  process  and/or  the  measurements  have  been  corrupted  by  additive  noise  (which 
is  a  random  process). 


Assuming  the  random  processes  are  all  normally  Gaussian  distributed,  and 
measurements  are  made  using  N  sensors,  the  associated  probability  density  function  is 
given  by  [5] 


0) 


1  -trace((X-M)Hc-x(x-M)) 

[Tr^detC]^ 


where  the  superscript  H  denotes  the  conjugate-transpose  operation,  and  the  vectors 
xo, x-k- l  represent  the  random  complex  measurement  data  associated  with  all  N 
sensors  for  time  instances  t  =  to,  h, ...,  tx-i  as  defined  by 


(2) 


x0(k) 
xi  (k) 
x2{k) 


for  0  <  k  <  K. 


XN-i(k) 


Additionally,  the  matrix  X  represents  all  K  measurement  vectors  as  given  by 


(3) 


X  =  [x0,xi,...,xtf_i], 


the  matrix  M  represents  the  corresponding  mean  values  of  the  measurements,  and  C  is 
the  N  x  N  covariance  matrix  describing  the  correlations  among  sensors.  The  exact 
definitions  of  the  matrices  M  and  C  depend  on  how  the  above  density  function  is 
applied  to  the  particular  estimation  problem  to  be  solved,  as  will  be  seen.  Once  these 
definitions  have  been  set  up,  the  maximum  likelihood  solution  is  found  by  maximizing 
the  probability  density  function  with  the  most  appropriate  valid  choices  of  M  and/or  C. 
These  maximizing  choices  for  M  and/or  C  can  then  be  related  back  to  the  signal 
parameters  of  interest. 

For  direction  finding  purposes,  xi , xk- i  are  taken  to  represent  the  complex 
baseband  outputs  from  an  array  of  N  antennas.  The  definitions  for  M  and  C  depend  on 
the  assumptions  made  about  the  signals.  Given  that  man-made  HF  signals  are  bounded 
(finite  power)  and  any  constant  modulus  properties  will  be  destroyed  by  the  dynamic 
nature  of  the  ionosphere  (e.g.  fading  and  Doppler  spreading),  it  is  reasonable  to  model 
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these  signals  as  being  zero-mean  stochastic  processes.  Hence  this  leads  to  the 
assumptions  that 

(4)  M  =  0 
where  0  is  a  matrix  of  all  zeros,  and 

(5)  Cmn  =  E[xm(k)xn(k)*]  for  0  <k<K  and  0<m,n<N 
Cmn  is  the  element  of  C  located  in  the  mth  row  and  the  nth  column. 


2.2  The  Cost  and  Model  Error  Functions 


The  probability  density  function  (1)  can  be  simplified  using  (4)  to  give 


(6) 


/(x0,x  = 


1 


[7rwdetC]^- 


0-trace(xHc-1x) 


The  objective  is  to  find  the  unknown  covariance  matrix  C  which  maximizes  this 
expression.  This  is  equivalent  to  maximizing  the  cost  function 
Lc  =  ln/(x0,xi,...,x*:_i)  or, 


(7)  Lc  =  —NK  ln(7r)  -  K  ln(det  C)  -  trace (X^C^X). 

Since  the  addition  or  multiplication  by  a  constant  value  has  no  effect  on  the 
maximization,  the  cost  function  can  be  simplified  to 

(8)  L  =  —  ln(det  C)  —  trace(RC_1) 
where  R  is  the  data  covariance  matrix  constructed  from  X  using 

(9)  R  =  ^-XXH. 


Although  not  required  for  the  theoretical  development  of  the  SPIRE  algorithm,  the  cost 
function  can  modified  to  provide  a  measure  of  the  modeling  error  which,  in  turn, 
provides  information  about  the  quality  of  the  estimates.  A  simple  model  error  function 
can  be  defined  as 

(10)  e  —  Lmax  —  L 

where  Lmax  is  the  cost  function  value  when  the  model  covariance  matrix  exactly 
matches  the  data  covariance  matrix,  or  C  =  R.  Expanding  the  expression  for  e  in  terms 
of  (8)  and  simplifying,  then 

(11)  e  ~  ln(det  C)  +  trace(RC-1)  —  ln(detR)  —  N 
The  usefulness  of  e  is  shown  later  on  in  Sections  4,  5  and  6. 
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2.3  Modeling  the  Signal  Environment 


The  modeling  aspect  comes  into  play  when  C  is  being  selected.  A  model  is  used  to 
generate  C  based  on  input  modeling  parameters  such  as,  for  example,  the  number  of 
signals,  signal  bearings,  signal  amplitudes/powers,  and  noise  powers.  For  this  reason, 

C  is  referred  to  as  the  model  covariance  matrix  in  this  report. 

The  procedure  for  determining  the  optimum  choice  for  C  (i.e.  the  choice  which 
maximizes  the  cost  function  L)  starts  by  choosing  initial  model  parameters,  generating 
the  corresponding  model  covariance  C,  and  then  determining  the  cost  function  value  L. 
The  model  parameters  are  then  successively  refined  and  C  recomputed  until  the 
maximum  cost  function  value  has  been  achieved.  The  model  parameter  values 
corresponding  to  this  maximum  value  are  then  taken  to  be  the  optimum  or  maximum 
likelihood  estimates. 


The  particular  model  used  to  generate  the  model  covariance  matrix  depends  on  many 
factors  including  the  transmitter(s)  and  receiver  characteristics,  the  signal  propagation 
environment,  the  noise  sources,  and  so  on.  One  way  to  set  this  model  up,  is  to  consider 
the  generation  of  synthetic  data  which  imitates  the  collected  data,  and  then  use  this 
synthetic  data  to  determine  C  according  to 

(12)  C  =  ^YYh 

where  Y  is  the  matrix  of  model  data  and  has  the  same  form  and  dimensions  as  X. 


Based  on  these  simplifying  assumptions,  the  received  signal  can  be  decomposed  as 
(13)  Y  =  Yi  +  Y2  +  ...  +  Ym  +  N 


where  the  matrices  Yi,  Y2, ...,  Y m  represent  the  model  data  for  the  M  individual 
signals,  and  N  is  the  modeled  noise.  The  matrix  Ym,  for  0  <  m  <  M,  can  be  defined 
in  vector  form  as 

(14)  Ym  =  eTOam 

where  em  is  the  steering  vector  (or  array  response  vector)  for  the  mth  signal,  and  aTO  is 
the  corresponding  signal  amplitude  vector.  The  definitions  for  the  elements  of  the 
steering  vector  are  given  by 


(15) 


g  j (10  sin  4>m  COS  7pm  +2/0  cos  (pm  cos  7pm  ) 
1  sin  4>m  cos  7pm  +yi  cos  <pm  cos  7pm ) 


gjyfls-l  sin /pm  COS  Ipm+VN-l  COS  (pm  COS  1pm) 


where  A  is  the  signal  wavelength,  xn  and  yn  are  the  Cartesian  coordinates  for  antenna  n 
(with  the  phase  center  of  the  array  located  at  the  origin),  4>m  is  the  azimuth  angle  of  the 
mth  signal  measured  clockwise  with  respect  to  the  Y-axis  of  the  coordinate  system,  and 
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tpm  is  the  elevation  angle  measured  with  respect  to  the  X-Y  plane  (the  ground).  The 
definition  for  the  elements  of  the  signal  amplitude  vector  is  given  by 


(16) 


<4(0) 

4(1) 


4(*-i) 


The  coefficients  om(0),  am(l), am(K  —  1)  are  the  received  complex  amplitudes  of 
the  mth  signal  for  time  instances  fo,  ti, tj<- i- 


The  model  covariance  can  also  be  written  as  a  sum  of  the  noise  covariances,  the  signal 
covariances,  and  the  signal  cross-covariances,  or 

M  M- 1  M 

(17)  C  =  cr2Cr?  +  ^  ]  C mm  +  ^  ^  ^  '  (Q mn  +  Q mn) 

m~  1  m—  1  n=m-f 1 

where  cr2  is  the  noise  power,  is  the  normalized  noise  covariance  matrix 
(trace  C,,  =  1),  Cmm  is  the  signal  covariance  matrix  for  the  mth  signal,  and  Qmn  is 
the  signal  cross-covariance  matrix.  The  generation  of  these  matrices  is  discussed  in  the 
following  paragraphs. 


The  noise  covariance  matrix  Cr;  is  assumed  to  be  known  a  priori  and  will  not  be 
considered  as  part  of  the  estimation  process.  The  determination  of  can  be  done 
either  through  theoretical  statistical  considerations,  or  through  measurements.  For 
example,  if  the  noise  is  known  to  be  white  Gaussian  in  nature  with  equal  but 
uncorrelated  amplitudes  in  each  channel,  then  E[j]m(k)r]:^l(k)]  =  cr2  and 
=  0  for  0  <  m,  n  <  N  and  m  ^n,  hence 

(18)  C„  =  jjIN 

where  I//  is  the  N  x  iV  identity  matrix.  Alternatively,  if  data  measurements  can  be 
taken  when  no  signals  are  present,  then 


(19) 


Cv  = 


XXH 
trace  (XX 


More  elaborate  procedures  could  be  developed  involving  a  number  of  measurement  sets 
with  the  same  noise  environment  but  different  signal  directions,  however  the 
development  of  this  kind  of  approach  is  beyond  the  scope  of  this  report.  It  suffices  to 
say  that  joint  estimation  of  both  the  noise  and  signal  characteristics  should  be  avoided  if 
possible  since  it  leads  to  poorer  accuracy. 


The  signal  covariance  matrix  Cmm  can  be  defined  in  terms  of  the  model  data  as 

(20)  cmm  =  -^Ymy£ 
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which  can  be  further  expanded  in  terms  of  the  component  amplitude  and  steering 
vectors  as 

nu  C  -  a™amc  eH 

In  a  similar  fashion,  the  cross-covariance  matrix  Qmn  can  also  be  defined  and 
expanded  to  get 

O  —  — Y 

Wmn  —  X  m  A  n 

(22)  -  -m^neme^  for  n  >  m. 

I\ 

For  the  high  latitude  HF  sky  wave  environment,  the  assumption  is  made  that  the 
cross-covariance  terms  disappear.  This  is  based  on  considering  the  various  cases. 
Signals  from  different  transmitters  will  be  uncorrelated.  Signals  originating  from  the 
same  transmitter  but  reflected  off  different  layers  of  the  ionosphere  will  also  be 
uncorrelated  due  to  the  large  path  length  differences  usually  encountered  (i.e.,  the  path 
delay  time  differences  will  be  greater  than  the  inverse  bandwidth  of  the  signal).  Signals 
originating  from  the  same  transmitter  but  scattered  from  different  parts  of  the  same 
region  (or  patch)  of  the  ionosphere  will  also  be  uncorrelated  since  the  scattering 
elements  within  the  region  are  short-lived  [11].  Given  Doppler  spread  measurements 
for  patch  scattered  signals  of  up  to  40  Hz  [10]),  it  is  further  assumed  that  this  is  caused 
by  the  birth/death  rate  of  the  scatterers  giving  a  decorrelation  time  of  r  >  l/40s.  More 
generally,  if  the  Doppler  spread  is  caused  by  turbulent  motion  of  the  electron  gas 
plasma  of  the  ionosphere,  from  the  perspective  of  the  receiver,  the  results  come  to  the 
same  thing  (i.e.  r  >  l/40s). 

A  violation  of  the  assumption  of  uncorrelated  signals  is  the  case  of  local  site  multipath. 
However,  like  the  problem  of  determining  the  noise  correlations,  including  the  signal 
correlations  in  the  estimation  process  is  highly  undesirable  (although  in  the  standard 
stochastic  ML  approach  this  is  done  [8]).  Hence  it  is  assumed  that  either  the  receiver 
site  is  well  chosen,  or  signal  correlations  can  be  determined  independently  and 
corrected  as  done,  for  example,  in  [9], 

The  end  result  is  that  (17)  can  be  simplified  to  become 

M 

(23)  c  =  cr^Cj)  +  y  ^  Cmm 

m~  1 

The  great  advantage  of  this  form  is  that  the  signal  part  of  the  model  covariance  matrix 
can  be  generated  as  the  sum  of  the  individual  signal  covariance  matrices.  More 
importantly,  by  dividing  the  regions  of  interest  into  sufficiently  small  subregions,  each 
subregion  can  be  represented  by  a  single  scattering  element  or  point  source  leading  to 

M 

(24)  C  =  <x2C„+J>m  eme^ 

m=  1 
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where  em  represents  the  steering  vector  associated  with  the  bearing  of  the  mth 
subregion  and  $m  is  the  corresponding  signal  power. 

The  optimum  solution  for  a2  and  sm  for  m  =  1, M  is  found  by  choosing  the  values 
which  maximize  the  cost  function  L.  If  the  region  to  be  subdivided  is  the  entire 
field-of-view  of  the  antenna  array,  then  the  result  will  be  a  maximum  likelihood 
estimate  of  the  spatial  power  spectrum  where  no  assumptions  have  been  made  about  the 
shape  and  power  distribution  of  the  scattering  region(s),  nor  any  assumption  about  the 
number  of  signals  (a  value  required  by  most  superresolution  methods). 

On  the  face  of  it,  given  M  will  be  very  large  (e.g.  for  a  1°  uniform  spacing  over  both 
0  —  360°  in  azimuth  and  0  —  90°  in  elevation  then  M  =  32400!),  there  would  appear  to 
be  two  fatal  objections  to  the  proposed  approach.  The  first  objection  is  that  the  model 
seems  to  be  over-determined  (too  many  model  parameters)  so  that  many  different  and 
incorrect  solutions  will  exist.  The  second  objection  is  that  estimation  of  so  many 
parameters  will  make  the  method  computationally  slow.  However,  the  fact  that  the 
model  parameters  are  constrained  so  that  a2,  sm  >  0,  plus  the  reduction  in  the  number 
of  kinds  of  parameters  (i.e.  power  only  versus  the  power,  bearing,  and  bearing  spread 
parameters  of  the  SML  algorithm)  has  a  major  impact  on  overcoming  these  objections 
as  will  be  seen  in  later  sections. 


DREO  TR  2000-099 


3.  MODEL  PARAMETER  ESTIMATION 


As  indicated  previously,  optimization  of  the  power  parameters  a 2  and  si, % 
proceeds  in  a  way  which  maximizes  the  cost  function  L.  The  basic  algorithm  begins  by 
initializing  the  noise  and  signal  power  parameters,  and  then  iteratively  updates  these 
parameters  using  a  gradient  technique  until  a  sufficiently  accurate  result  is  achieved. 
The  basic  algorithm  is  listed  below.  Additional  details  for  some  of  the  steps  listed  can 
be  found  in  the  indicated  sections. 

1 .  Initialize  model  noise  power  by  performing  an  eigendecomposition  on  the 
whitened  data  covariance  matrix  to  get 

_a  if 

(25)  <VR<V 

i—  1 

and  then  overestimating  the  initial  value  of  cr2  using 

(26)  cr2  =  2Ajv 
(Section  3.1). 

2.  Initialize  the  model  signal  powers  by  setting 

(27)  Si  =  S2  =  •  •  •  =  SM  =  -^(traceR  -  a2) 

(Section  3.2). 

3.  Initialize  the  loop  counter:  loop  =  0 

4.  Increment  the  loop  counter:  loop  — »  loop  +  1 

5.  Update  signal  power  parameter  estimates  (Section  3.3). 

6.  While  loop  <  maxJoop  go  to  step  4. 

7.  (Optional)  Update  noise  power  parameter  (Section  3.4). 

8.  Output  model  parameter  estimates. 


A  value  of  maxJloop  =  20  has  been  found  to  give  good  results.  Enhancements  to 
accelerate  the  processing  speed  and  appropriate  choices  for  the  various  algorithm 
control  parameters  are  discussed  in  Section  4. 
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3.1  Initial  Noise  Power  Estimate 


The  initial  estimate  of  the  noise  power  a2  is  based  on  the  idea  that,  in  the  presence  of 
spatially  white  Gaussian  noise,  the  data  covariance  matrix  can  be  divided  into  signal 
and  noise  subspaces.  In  the  derivations  thus  far,  the  noise  has  not  been  assumed  to  be 
white,  but  it  has  been  assumed  to  be  known.  Hence  the  first  step  is  to  whiten  the  noise 
by  performing  the  following  modification  to  the  data  covariance  matrix. 

(28)  Rw  =  WhRW 
where 

(29)  WWH  =  C"1 

There  are  many  equally  valid  solutions  for  W  which  satisfy  (29)  but  for  simplicity  of 
notation  the  choice 

(30)  W  = 

has  been  used.  The  actual  choice  is  a  matter  of  convenience.  Additionally,  in  the  special 
case  where  the  noise  is  already  white  and  Gaussian,  then 

(31)  Rw  =  NR 


Using  eigendecomposition,  the  whitened  data  covariance  matrix  can  be  represented  by 

N 

(32)  Rw  =  ^2  AiVjvf 

i= 1 

where  Ai, ...,  Ay  are  the  eigenvalues  ordered  so  that  Aj  >  •  •  •  >  Ay,  and  vi,  ...vy  are 
the  corresponding  orthonormal  eigenvectors.  In  the  ideal  case,  R„  is  formed  from  an 
infinite  number  of  sensor  snapshots,  there  are  M  <  N  point-source  signals  impinging 
on  the  array,  and  the  measurements  are  corrupted  by  additive  white  Gaussian.  Under 
these  conditions,  the  first  M  eigenvectors  will  be  associated  with  signal  +  noise,  while 
the  rest  (N  —  M  eigenvectors)  will  be  associated  with  noise  only.  The  corresponding 
values  of  the  N  —  M  smallest  eigenvalues  will  all  equal  a2.  Under  these  ideal 
conditions,  estimating  cr2  from  any  of  the  smallest  N  —  M  eigenvalues  is  a  trivial 
exercise. 

In  the  high  latitude  HF  case,  the  number  of  measurements  will  not  be  infinite,  and  the 
signals  will  not  be  point-source  (due  to  azimuth/elevation  spreading).  For  a  limited 
number  of  point-source  signals,  the  problem  of  finite  samples  could  be  solved  by 
averaging  the  smallest  N  —  M  eigenvalues  to  produce  an  estimate  of  the  noise  power. 
The  problem  of  spatial  spreading  cannot  be  solved  so  easily,  however,  since  spreading 
causes  the  effective  value  of  M  (remembering  that  M  is  the  number  of  point-source 
signals)  to  increase  so  that  it  may  exceed  N.  Figure  1  illustrates  this  point  through 
simulation.  Figure  la  shows  the  spatial  power  spectrum  consisting  of  a  point-source 
signal,  a  spread-source  signal,  and  a  noise  floor  set  to  yield  a  signal-to-noise  ratio 
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(SNR)  of  20  dB.  Figure  lb  shows  the  corresponding  eigenvalues  for  three  cases:  signal 
only;  noise  only;  and  signal  plus  noise.  The  antenna  array  configuration  used  for  the 
simulation  is  shown  in  Figure  2  and  discussed  in  Section  4. 1 .  The  results  demonstrate 
that,  in  this  case,  the  signal  power  affects  all  the  eigenvalues  implying  M  >  N  even 
thought  there  were  only  two  signals. 

The  simplest  solution  to  the  spreading  problem  is  to  assume  that  at  least  one  of  the 
dimensions  of  Rw  is  dominated  by  noise.  The  noise  power  can  then  be  estimated  as, 

(33)  <r2  =  XN 

For  the  example  in  Figure  1,  choosing  the  smallest  eigenvalue  for  the  noise  power 
estimate  gives  the  best  result  even  though  it  is  slightly  overestimated. 

In  cases  where  the  signal  power  dominates  all  dimensions  of  R^  (e.g.,  in  the  previous 
example  this  would  be  true  for  an  SNR  =  30  dB),  then  \N  will  be  a  large  overestimate 
of  a2.  In  practice,  this  overestimate  has  not  been  found  to  significantly  affect  the 
accuracy  of  the  signal  estimates  (which  are  of  most  interest),  so  no  attempt  has  been 
made  to  investigate  this  problem  further. 

In  fact,  it  has  been  found  that  purposely  overestimating  the  noise  power  by  a  factor  of 
two,  and  then  fixing  this  value  until  after  the  final  iteration  of  the  SPIRE  algorithm, 
leads  to  more  accurate  estimates  of  the  noise  power  and  better  detection  of  weaker 
signals.  If  the  noise  power  is  updated  within  the  loop,  the  signal  model  sometimes 
adapts  to  the  noise  (due  to  the  signal  model’s  greater  flexibility)  forcing  the  estimated 
noise  power  to  zero  and  causing  false  signal  estimates.  Under  these  circumstances 
convergence  occurs  on  an  undesirable  false  maximum  of  the  cost  function.  Hence  it  is 
better  to  wait  until  the  signal  model  estimate  has  stabilized  before  fine  tuning  the  noise 
power. 


3.2  Initial  Signal  Power  Estimate 

The  initial  signal  powers  are  set  equal  according  to 

(34)  si  =  «2  =  •  •  •  =  sm  =  s 

Since  the  total  power  observed  in  the  data  is  assumed  to  be  produced  by  uncorrelated 
processes,  then  ideally 

M 

(35)  traceR  =  a2  +  sm 

m= X 

Hence,  given  the  signal  powers  are  all  equal,  the  quantity  s  is  given  by 

(36)  s  =  if  (^ceR  ~  ^ 

where  the  value  for  cr2  used  here  is  the  initial  value  computed  according  to  (33). 
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Eigenvalue  Number 

(b) 


Figure  1:  Simulation  of  a  spread-source  signal  environment  showing  (a)  the  spatial  power  spectrum,  and  (b)  the 
corresponding  eigenvalues  when  an  array  with  12  antennas  is  used. 
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For  a  limited  sample  of  sensor  data,  (36)  will  only  be  approximately  true  and  a  more 
accurate  value  of  s  could  be  obtained  by  maximizing  the  cost  function  L.  However,  for 
the  sake  of  initialization,  and  in  the  interests  of  reducing  the  number  of  computations, 

(36)  is  sufficient. 

3.3  Signal  Power  Estimates 

The  optimum  solution  for  the  signal  power  sm  can  be  found  by  setting  the  gradient  of  L 
to  zero  where  the  gradient  is  measured  with  respect  to  sm.  The  generic  gradient  of  the 
cost  function  L  (derived  in  [3])  is  given  by 

(37)  G(a)  =  trace  ((Cr'R  - 

In  the  case  where  a  =  sm,  setting  the  gradient  G(sm)  to  zero  leads  to  the  following 
expression  (based  on  the  derivation  in  [3]  with  some  slight  rearranging), 

e^(C-1R-I)C-1em 

(38)  m  ~  (e£C-iem)2  ' 

where  A sm  is  the  estimated  signal  update.  Hence  for  a  particular  signal  source  m, 
sm  +  A sm  is  the  updated  signal  power  which  maximizes  L. 

For  computational  efficiency,  the  signal  power  updates  for  m  =  1, ...,  M  are  done  in 
parallel  so  that  the  matrices  CT1  and  (C-1R  -  I)C_1  need  only  be  calculated  once, 
and  not  M  times  for  each  individual  signal  power  update.  However,  this  results  in  an 
overestimation  of  the  signal  updates  so  a  compensating  scaling  factor  is  required. 
Calling  the  scaling  factor  //,  the  update  is  applied  as.  ^ 

To  ensure  the  fastest  convergence  using  this  approach,  the  value  of  n  which  maximizes 
L  is  calculated  for  each  iteration.  As  no  direct  solution  has  been  found,  the  bisection 
method  has  been  chosen  as  providing  a  reasonably  fast  search  procedure  to  determine 
the  desired  solution.  Using  the  generic  search  parameter  a,  this  method  is  shown  below. 

1.  Set  the  upper  and  lower  search  bounds,  aiow  and  a  high- 

2.  Estimate  the  search  parameter  as:  a  =  \{otiow  +  &high ) 

3.  Compute  the  gradient  (7(0:) . 

4.  If  G(a)  >0  then  atow  =  a,  otherwise  ahigh  =  a. 

5.  Repeat  steps  2-5  nine  more  times. 

Steps  2-5  can  be  repeated  more  times  for  greater  accuracy,  however  in  practice,  the 
number  of  times  indicated  has  been  found  to  be  sufficient  for  the  use  of  the  SPIRE 
algorithm. 
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The  appropriate  model  for  the  updated  covariance  matrix  is 

M  M 

(39)  C  =  a2C„+5>m  eme^  +  m  ^  Asmeme^ 

m—\  m— 1 

Note  that  the  matrix  formed  from  the  first  two  terms  on  the  right  hand  side  of  this 
expression  represents  the  model  covariance  before  updating  and  only  needs  to  be 
calculated  once,  not  in  every  iteration  of  the  bisection  method.  Using  this  model  and 
solving  (37)  with  respect  to  //,  the  corresponding  gradient  function  is  given  by 

M 

(40)  G(jx)  =  ]T  Asme^(C-1R  - 

i=m 

Finally,  the  search  limits  for  /i  are  given  by:  plow  =  0  and  n high  —  1. 

In  the  event  that  sm  +  fJ.Asm  <  0,  sm  is  assigned  a  value  of  zero  since  the  signal 
powers  can  never  be  negative.  Additionally,  in  the  event  that  s*  =  0  and  As*  <  0,  the 
calculation  involving  m  =  k  in  (39)  can  be  ignored  reducing  the  overall  number  of 
computations. 

There  are  two  problems  with  the  bisection  method  as  described  here:  the  function  G(ji) 
is  not  guaranteed  to  be  monotonically  decreasing  in  the  region  of  interest;  and  the 
method  does  not  take  into  account  the  constraint  sm  +  /j,A sm  >  0.  In  practice, 
provided  that  calculations  associated  with  obviously  negative  powers  are  eliminated 
(i.e.  Sk  =  0  and  Ask  <  0),  the  bisection  method  returns  an  answer  for  )i  which, 
although  not  always  optimum,  is  sufficiently  good  for  practical  purposes. 

3.4  Noise  Power  Update 

The  noise  power  update  a 2  can  also  be  calculated  using  the  bisection  method  outlined 
in  Section  3.2.  The  model  covariance  defined  in  (24)  is  suitable  for  these  purposes. 

This  definition  is  repeated  here  as 

M 

(41)  C  =  a2Cv  +  sm 

2—1 

Using  this  model  and  solving  (37)  with  respect  to  a2,  the  gradient  function  is  given  by 

(42)  G{a2)  =  trace  ((C_1R  -  I)C_1C^) 

The  lower  search  limit  for  a 2  is  afow  =  0.  The  upper  search  limit  is  given  by 

(43)  a\igh  -  trace  R 

where  trace  R  is  an  estimate  of  the  signal  plus  noise  power  (and  hence  an  overestimate 
for  the  noise  power). 

As  discussed  earlier  in  Section  3. 1 ,  the  noise  power  update  is  only  applied  after  the 
signal  model  estimates  have  been  refined  and  only  if  a  more  accurate  estimate  of  the 
noise  power  is  desired. 
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4.  IMPLEMENTATION  ISSUES 


To  investigate  some  of  the  practical  issues  associated  with  the  implementation  of  the 
SPIRE  algorithm,  a  series  of  simulations  were  carried  out.  The  details  of  the 
simulations  are  discussed  in  the  next  section  (Section  4.1).  The  implementation  issues 
addressed,  including  spatial  ambiguities,  signal  model  grid  spacing,  data  sample  size, 
and  limitation  on  the  number  of  signals,  are  discussed  in  Sections  4.2-4.5. 


4.1  Simulations  Settings 

For  these  simulations,  the  antenna  array  shown  in  Figure  2  was  chosen,  as  this  array 
geometry  was  investigated  in  [16],  [4]  and  was  found  to  have  very  good  characteristics 
for  direction  finding.  Assuming  an  ideal  free  space  response,  then  for  any  given  signal 
bearing,  the  azimuth  beamwidth  of  the  array  is  relatively  constant  with  respect  to  the 
azimuth  bearing  but  varies  with  the  elevation  bearing  according  to 

(44)  4>bw  =  7.8°/|cos,0|  for  \xp\  <  90° 

where  the  beamwidth  is  defined  as  the  angular  width  of  the  main  lobe  of  the  antenna 
array  gain  pattern  measured  at  the  3  dB  points  (i.e.  the  points  at  0.707 x  the  maximum 
gain).  Note  that  measuring  the  azimuth  beamwidth  at  or  near  xp  =  90°  is  effectively 
meaningless. 

The  elevation  beamwidth  is  given  by 

(45)  ipBW  —  7.8°/|  sin^l  for  \xp\  >  30°. 

For  elevation  bearings  below  30°,  the  beamwidth  is  somewhat  more  complicated  as 
shown  in  Figure  3.  The  failure  of  the  above  expression  at  the  lower  elevation  angles  is 
due  to  the  distortion  of  the  main  lobe  in  the  antenna  pattern  at  the  elevation  angle  xp  by 
the  reflection  of  this  lobe  at  the  elevation  angle  —ip.  For  example,  a  2-dimensional  x-y 
array,  using  the  free  space  assumptions,  has  a  symmetrical  gain  pattern  for  elevation 
angles  above  and  below  the  horizon,  i.e.,  a  main  lobe  and  a  reflection  lobe.  At  low 
signal  elevations,  the  main  and  reflection  lobes  begin  to  join.  They  are  considered 
merged  when  the  minimum  gain  between  the  two  lobes  is  greater  than  0.707 x  the 
maximum  gain.  Using  the  array  configuration  shown,  this  occurs  at  21.2°.  At  this 
point,  the  beamwidth  effectively  doubles.  For  even  lower  signal  elevations  the  merged 
lobes  move  closer  together  so  the  beamwidth  actually  decreases. 

The  simulation  signal  model  was  set  up  using  two  signals  with  the  spatial  power 
spectrum  shown  in  Figure  4.  The  spread-source  at  {(p,xp)  =  (10°,  30°)  had  a  total 
power  three  times  that  of  the  point-source.  The  noise  power  level  (represented  by  the 
floor  at  0  dB  in  Figure  4)  was  set  to  be  20  dB  less  than  the  point-source  signal  power. 

The  basic  steps  of  the  SPIRE  algorithm  were  outlined  in  the  previous  section. 
Following  these  steps,  the  result  of  processing  simulated  data  using  a  subregion  size  of 
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Figure  2:  Three  dimensional  view  of  the  antenna  array.  Each  grid  square  has  a  dimension  of  1 A  x  1  A. 


Figure  3:  Antenna  array  elevation  beamwidth  as  a  function  of  elevation  comparing  the  simulated  response  (solid 
line)  with  the  sin- 1  ip  predicted  response  (dashed  line). 
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Figure  4:  Spatial  power  spectrum  of  signal  model  used  for  simulations. 


1°  x  1°  (i.e.  the  model  covariance  matrix  was  constructed  from  point  source  signals 
placed  every  degree  in  azimuth  and  every  degree  in  elevation)  and  a  sample  size  of 
K  =  10000  is  shown  plotted  in  Figure  5a.  Despite  the  large  amount  of  data  used  to 
generate  the  results,  the  shape  of  the  SPIRE  estimate  power  spectrum  is  less  than  ideal. 
Improvements  to  the  shape  are  discussed  in  Section  4.2. 

For  comparison  purposes,  the  Minimum  Variance  (MV)  [12]  and  MUSIC  [14] 
algorithms  were  also  used  to  process  the  results  and  the  results  are  plotted  in  Figure  5b 
and  c.  A  search  grid  employing  the  same  1°  x  1°  spacing  as  used  for  the  SPIRE 
algorithm,  was  also  used  by  both  the  MV  and  MUSIC  algorithms.  Since  MUSIC  was 
originally  developed  for  direction-of-arrival  estimation,  and  not  spectrum  estimation, 
the  square  root  of  the  MUSIC  output  has  been  displayed  here  and  all  other  MUSIC  plots 
in  this  report  since  this  results  in  a  better  estimate  of  the  relative  spectral  power  levels. 

The  processing  time  for  the  example  shown  in  Figure  5  was  38  seconds  for  basic 
SPIRE  and  1.6  seconds  each  for  MV  and  MUSIC.  A  faster  method  of  producing  the 
SPIRE  results  is  discussed  in  the  section  4.3. 


4.2  Spatial  Ambiguities 

Given  the  large  amount  of  data  used  to  generate  the  data  covariance  matrix 

(K  =  10000),  a  better  correspondence  might  have  been  expected  between  the  actual 
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signal  model  shown  in  Figure  4  and  the  SPIRE  estimated  model  shown  in  Figure  5a. 
However,  the  results  illustrate  an  apparent  ambiguity  problem  with  the  spatial  power 
spectrum.  For  example,  nearly  identical  results  to  those  shown  in  Figure  5  can  be 
produced  with  the  signal  model  having  the  spatial  power  spectrum  shown  in  Figure  6. 
In  other  words,  the  two  signal  models  produce  almost  identical  data  covariance 
matrices  although  the  spread-source  signal  regions  have  different  spatial  shapes. 
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Figure  6:  Spatial  power  spectrum  of  alternate  signal  model  which  leads  to  virtually  the  same  results  as  in  Figure  5 
when  the  model  generated  data  is  processed  with  basic  SPIRE,  MV,  or  MUSIC. 

This  result  is  not  totally  unexpected  since  as  smaller  and  smaller  features  are 
considered,  the  angular  resolution  limitations  of  the  antenna  array  will  become  a 
problem.  Hence,  an  arbitrarily  shaped  region  will  be  indistinguishable  from  a  point 
source,  or  any  other  arbitrarily  shaped  region  for  that  matter,  as  long  as  the  regions  are 
very  small  and  at  the  same  angular  position.  Applying  this  idea  to  the  features  in 
Figures  4  and  6,  “small”  means  less  than  the  antenna  array  beamwidth. 

Given  this  resolution  problem,  it  makes  more  sense  to  estimate  regions  using  the 
simplest  shape  that  yields  an  acceptable  solution  rather  than  more  complex  shapes. 
Using  this  as  a  guideline,  and  investigating  a  number  of  different  modifications  to  the 
basic  SPIRE  algorithm,  the  following  modified  algorithm  was  developed. 

1.  Perform  basic  SPIRE  as  outlined  in  Section  3  except  using  maxJoop  —  10. 

2.  Subdivide  the  signal  grid  into  P  regions. 
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3.  Do  the  following  steps  4-12  twice. 

4.  Initialize  the  region  counter:  p  =  0. 

5.  Increment  the  region  counter:  p  — t  p  +  1. 

6.  Select  only  grid  signals  from  region  p  for  further  processing. 

7.  Reduce  all  signal  powers  in  region  p  by  a  factor  of  2. 

8.  Do  steps  9-11  maxJtoop  times. 

9.  Include  grid  signals  bordering  region  p. 

10.  Update  region  m  signal  power  estimates  (Section  3.3). 

1 1  •  Remove  any  grid  signals  with  zero  power. 

12.  While  p  <  P  go  to  step  5. 


After  running  the  basic  SPIRE  algorithm,  the  signal  grid  will  consist  of  mainly  zero 
power  signals  with  islands  or  regions  of  positive  power  signals.  The  first  step  in  the 
modification,  then,  is  to  identify  these  regions.  Once  this  is  done,  each  region  is 
separately  processed  using  the  basic  SPIRE  approach  except  that  signal  powers  are 
initialized  to  half  their  previous  value,  and  the  regions  are  allowed  to  grow  and  shrink  in 
size  (steps  9  and  11). 

Subdividing  the  grid  into  regions  reduces  the  amount  of  processing  since  the  “zero” 
parts  of  the  signal  grid  are  ignored.  It  also  improves  convergence.  Halving  the  signal 
powers  results  in  simpler  spatial  shapes.  Steps  9  and  1 1  also  allow  the  region  to  develop 
in  a  smooth  manner. 

Using  this  enhanced  form  of  SPIRE  for  the  same  data  used  to  produce  Figure  5,  the 
result  is  shown  in  Figure  7.  Although  not  identical  to  Figure  4,  the  results  are 
considerably  improved  compared  to  Figure  5  when  the  “simpler  is  better”  guideline  is 
followed. 

An  upper  limit  may  also  be  imposed  on  the  number  of  regions  P  to  further  reduce 
processing  without  significantly  affecting  the  results  as  long  as  P  equals  or  exceeds  the 
true  number  of  spread-source  signals. 


4.3  Signal  Model  Grid  Spacing 

The  angular  spacing  between  point-source  signals  used  in  the  signal  model  grid  must  be 
sufficiently  narrow  for  the  model  to  represent  a  real  spread-source  signal  adequately. 
However,  making  this  spacing  too  narrow  can  unnecessarily  increase  the  number  of 
computations  since  the  increase  in  computations  is  inversely  proportional  to  the  square 
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Figure  7:  Results  of  processing  K  =  10000  samples  of  simulated  data  using  the  enhanced  form  of  SPIRE. 


of  the  spacing.  Hence  it  is  useful  to  determine  the  widest  spacing  that  can  be  used 
before  introducing  too  much  error. 

Since  the  resolving  power  of  any  antenna  array  is  a  function  of  the  beamwidth  of  the 
array,  a  natural  choice  would  be  to  use  a  spacing  value  which  is  some  fraction  of  the 
minimum  beamwidth  of  the  array.  Figure  8  shows  several  examples  of  the  results  of 
processing  the  same  simulation  example  as  before  (see  Figure  5f)  where  the  signal 
model  grid  spacing  was  varied  from  20%  to  100%  of  the  beamwidth  (1.6°  to  7.8°). 

Comparing  the  different  spacings,  good  results  were  obtained  for  spacings  up  to  80%  of 
the  beamwidth.  At  80%  or  more,  extraneous  signals  begin  to  appear  in  the  results 
suggesting  that  this  represents  the  upper  limit  on  the  spacing.  A  similar  conclusion  was 
also  drawn  in  [3]  for  the  SML  algorithm  based  on  an  analysis  of  bearing  accuracy 
versus  grid  spacing. 

An  advantage  of  the  larger  spacings  is  that  since  the  processing  time  is  inversely 
proportional  to  the  square  of  the  spacing,  larger  spacings  means  faster  processing.  A 
problem  with  larger  spacings,  however,  is  illustrated  in  Figure  9.  In  this  example,  two 
closely  spaced  point-source  signals  are  detected  but  not  resolved  for  spacings  that  are 
greater  than  half  the  angular  distance  between  the  signals  (in  this  case  40%  of  the 
beamwidth  or  more).  Hence  one  way  to  reduce  the  amount  of  processing,  yet  maintain 
acceptable  resolution,  would  be  to  employ  larger  grid  spacings  for  the  early  iterations 
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Figure  8:  Results  for  SPIRE  using  various  signal  grid  spacings  showing  (a)  the  simulation  signal  model,  and  the 
estimated  models  for  spacings  of  (b)  20%,  (c)  40%,  (d)  60%,  (e)  80%,  and  (f)  100%  of  the  antenna  array 
beamwidth. 
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and  then  finer  spacings  for  the  later  iterations. 


(c) 


(d) 


Figure  9:  Resolution  of  two  point-source  signals  with  bearings  (<j> 1,^1)  =  (-1.5°, 29°)  and 

((p2,ip2)  =  (1.5°, 26°)  and  signal  powers  sj/a 2  =  20  dBands\jc r2  =  30  dB.  The  SPIRE  signal  grid  spacings 

used  were  (a)  1°,  (b)  2°,  (c)  3°,  and  (d)  4°. 

Applying  this  coarse/fine  spacing  idea  to  the  SPIRE  algorithm,  a  fast  version  has  been 
developed.  Using  the  spacing  scheme  illustrated  in  Figure  10,  the  first  10  iterations  are 
carried  out  using  the  coarse  grid  in  the  same  manner  as  before  but  without  the  final 
noise  power  update.  The  signal  power  estimates  are  then  interpolated  from  the  coarse 
grid  to  the  fine  grid  using  a  simple  quadratic  function.  For  example,  all  the  signal 
powers  of  the  coarse  grid  shown  in  the  dotted  box  in  Figure  10  are  used  to  interpolate 
the  fine  grid  signal  powers  within  the  solid  box  based  on  the  expression 

(46)  s  =  co  +  c\x  +  c2y  +  c^xy  +  c±x2  +  c5y2 

where  x,y  =  {  —  1,0, 1}  are  the  Cartesian  coordinates  of  the  signals  in  the  solid  box 
(the  spacing  between  adjacent  fine  grid  signals  is  1  unit  and  the  center  of  the  box 
represents  the  origin).  The  coefficients  cq,  ...,05  are  determined  using  least  squares 
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estimation  techniques  in  conjunction  with  the  known  coarse  grid  signal  powers  and 
x,  y  —  {—3, 0, 3}.  The  coefficients  are  then  divided  by  9  to  maintain  the  same  signal 
power  density  since  the  fine  grid  contains  9  times  as  many  model  signals.  The 
coefficients  are  further  divided  by  a  factor  of  2  which  causes  an  underestimate  of  s  in 
(46),  but  ultimately  results  in  a  smoother  transition  from  the  coarse  to  the  fine  grid 
signal  power  estimate. 


Azimuth 

Figure  10:  Grid  layout  for  coarse/fine  processing  scheme.  The  large  dots  represent  the  positions  of  the  point-source 
signals  used  for  initial  (coarse)  processing  and  the  small  circles  plus  large  dots  represent  the  positions  used  for  final 
(fine)  processing.  The  coarse  grid  signals  within  the  dotted  box  are  used  to  interpolate  the  fine  grid  values  within  the 
solid  box. 

Using  the  pared  fine  grid,  ten  more  iterations  are  then  carried  out  followed  by  the  final 
update  of  the  noise  power. 

Along  with  the  coarse/fine  processing  scheme,  a  further  reduction  in  the  amount  of 
processing  can  be  achieved  by  only  updating  signals  in  the  fine  grid  which  either  have 
positive  (nonzero)  powers  or  border  signals  that  have  positive  power.  Since  signals  will 
usually  be  limited  to  small  regions  of  the  sky,  the  reduction  in  processing  can  be 
substantial.  For  example,  using  a  3°  coarse  spacing  and  1°  coarse  spacing  for  the  signal 
model  grid,  the  computation  time  needed  to  produce  Figure  5a  was  reduced  by  a  factor 
of  almost  8  (5.0  seconds  versus  37.9  seconds). 

Using  the  coarse/fine  scheme  for  the  enhanced  form  of  SPIRE  discussed  in  Section  4.2, 
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instead  of  the  basic  form,  producing  Figure  5f  resulted  in  a  reduction  in  the 
computation  time  of  a  factor  of  slightly  more  than  3  (1 1.2  seconds  versus  37.9 
seconds).  The  smaller  improvement  in  this  case  was  due  to  the  fact  that  a  portion  of  the 
computation  time  was  taken  up  by  extra  processing  stage  which  is  the  same  whether  the 
fast  or  basic  version  of  SPIRE  is  used. 

No  doubt  variable  spacing  schemes  could  be  successfully  employed  to  reduce 
processing  further  (e.g.,  spacing  as  a  function  of  beamwidth  and/or  spherical  distances 
between  grid  points),  but  this  avenue  was  not  investigated.  Additionally,  since  the 
computations  of  the  signal  powers  can  be  performed  in  parallel,  much  faster  versions  of 
the  SPIRE  algorithms  could  be  implemented  using  parallel  processors. 


4.4  Data  Sample  Size 

As  stated  previously,  the  key  assumption  in  the  development  of  the  SPIRE  algorithm  is 
that  all  signals  are  uncorrelated  with  each  other  including  signals  from  different  parts  of 
the  same  spread  region.  Consequently  for  collection  purposes,  a  sufficient  number  of 
data  samples  must  be  collected  to  ensure  that  proper  decorrelation  occurs  in  the  sample 
data. 

To  investigate  the  effect  of  the  number  of  data  samples  K,  a  series  of  simulations  was 
carried  out  using  the  same  signal  model  used  previously  (shown  in  Figure  4),  but 
varying  the  number  of  samples.  The  results  are  shown  in  Figure  1 1 . 

The  results  show  that  the  spread  region  is  most  affected,  becoming  “hillier”,  as  the 
number  of  samples  decrease.  This  is  a  result  of  the  sensitivity  of  the  estimation  process 
to  the  kinds  of  ambiguities  discussed  in  Section  4.2.  Small  perturbations  due  to  noise 
and  incomplete  decorrelation  can  have  large  effects  on  the  estimated  spatial  shape  of 
the  spatial  region.  Since  larger  spread  regions  require  a  greater  number  of  grid  signals 
for  modeling,  there  is  a  greater  possibility  for  ambiguities.  By  way  of  comparison,  a 
point-source  signal  only  requires  one  grid  signal  with  no  possibility  of  any  ambiguity. 
For  this  particular  set  of  results,  the  spread  region  is  accurately  modeled  for  K  >  1000 
while  the  point-source  signal  is  accurately  modeled  for  all  values  of  K,  although  the 
estimated  signal  exhibits  spreading  in  elevation  at  the  smaller  values. 

A  second  series  of  simulations  was  carried  out  featuring  the  “H”  shaped  signal  model 
shown  in  Figure  12.  As  before,  the  number  of  samples  was  varied.  The  effect  of 
decreasing  the  number  of  samples  in  this  case  is  similar  to  the  previous  example. 

To  illustrate  the  improved  abilities  of  the  SPIRE  algorithm  to  estimate  the  spatial  power 
spectrum  versus  MV  or  MUSIC,  some  comparative  results  are  shown  in  Figure  13  for 
K  =  1000.  The  MV  and  MUSIC  results  shown  in  Figure  13  e  and  f  were  enhanced  by 
showing  only  the  spectrum  within  10  dB  of  the  peak  value.  This  10  dB  threshold  was 
based  on  knowledge  of  the  actual  spatial  power  spectrum  of  the  signal  -  knowledge 
which  wouldn’t  normally  be  available  in  practical  applications. 
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Figure  1 1:  The  effect  of  varying  the  number  of  data  samples  on  the  performance  of  SPIRE.  The  results  shown  are 
for  (a)  K  =  10000,  (b)  K  =  1000,  (c)  K  =  300,  (d)  K  =  100,  (e)  K  =  50,  and  (f)  K  =  25. 


26 


DREO  TR  2000-099 


CD 

"D 


CD 

T3 


(e) 


(0 


Figure  13:  A  comparison  of  the  estimation  accuracy  of  SPIRE,  MV,  and  MUSIC  for  the  “H”  shaped  signal  using 
K  =  1000  and  showing  (a)  the  true  signal  spatial  power  spectrum,  (b)  the  SPIRE  result,  (c)  the  MV  result,  (d)  the 
MUSIC  result,  (e)  the  modified  MV  result,  and  (f)  the  modified  MUSIC  result  For  the  modified  results,  only  features 
within  10  dB  of  the  maximum  peak  value  are  shown . 
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Finally,  the  effect  of  sample  size  on  the  measured  model  error  e  is  shown  in  Figure  14. 
The  plotted  model  error  is  the  average  of  the  error  values  generated  when  processing 
the  simulated  data  used  to  produce  Figure  1 1  (the  spread-source  and  point-source 
signals)  and  Figure  12  (the  “H”-shaped  signal).  The  results  show  that  the  “H”-shaped 
signal  was  easier  to  model  than  the  spread-source  signal  plus  point-source  signal.  The 
reasons  are  likely  due  to  the  ambiguity  problem  discussed  earlier  (the  wider  angular 
area  of  the  spread-source  signal  would  suffer  the  ambiguity  problem  more  than  the 
narrow  ridges  of  the  “H”-shaped  signal),  and  the  quantization  effect  (the  point-source 
can  only  be  exactly  modeled  if  its  bearing  corresponds  with  the  bearing  of  a  signal  in 
the  signal  model  grid). 


Figure  14:  Model  error  as  a  function  of  sample  size.  The  upper  curve  (solid  line)  shows  the  results  for  the  same  data 
used  in  Figure  1 1  (point-source  and  spread-source  signals).  The  lower  curve  (dashed  line)  shows  the  results  for  the 
same  data  used  in  Figure  12  (“H”-shaped  signal). 

4.5  Number  of  Signals 

The  general  rule  used  to  determine  the  maximum  number  of  signals  that  may  be 
estimated  using  a  superresolution  algorithm  is  N  —  1  where  N  is  the  number  of 
sensors.  This  rule  applies  to  algorithms  which  make  no  assumption  about  signal 
correlation  (i.e.  the  correlations  are  implicitly  or  explicitly  estimated). 

The  SPIRE  algorithm  assumes  that  the  signals  are  fully  uncorrelated  so  that  estimation 
of  the  signal  correlations  is  not  required.  As  a  result,  the  N  —  1  limit  can  be  exceeded 
as  illustrated  in  Figure  15  where  the  bearings  of  13  uncorrelated  point-source  signals 
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(Figurel5a  and  b)  are  correctly  estimated  using  an  array  of  12  antennas  (Figurel5c  and 
d).  A  few  spurious  signals  were  also  estimated  with  the  maximum  false  peak  in  the 
spectrum  still  5  dB  below  the  minimum  true  signal  peak. 

The  MV  algorithm  also  has  the  capability  to  exceed  the  N  —  1  limit  as  illustrated  in 
Figure  15e  and  f  where  the  13  largest  peaks  correspond  to  the  true  signals.  However,  the 
resultant  spectrum  is  considerably  more  difficult  to  interpret  than  the  SPIRE  spectrum 
with  many  more  false  peaks,  and  with  some  of  them  coming  within  1  dB  of  a  true  peak. 

The  MUSIC  algorithm  is  limited  to  N  —  1  signals  and  consequently  no  results  are 
shown  here. 

Although  the  results  suggest  that  there  is  no  hard  upper  limit  for  uncorrelated  signals, 
they  illustrate  that  estimating  an  increasing  number  of  signal  directions  comes  at  the 
price  of  lower  accuracy.  The  SPIRE  algorithm  was  able  to  estimate  all  13  signals  in  the 
example,  but  the  results  are  not  as  good  as  the  results  shown  in  Figure  9a  or  Figure  1  lb 
where  only  two  signals  were  involved.  From  the  practical  point-of-view,  it  seems  more 
reasonable  to  use  the  SPIRE  algorithm  for  environments  where  the  number  of  signals  is 
somewhat  less  than  the  number  of  sensors. 
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Figure  15:  An  example  of  estimation  of  13  point-source  signals  using  12  antennas,  K  =  1000,  and  SNR  =  20  dB, 
showing  (a)  the  actual  spatial  signal  spectrum,  and  (b)  azimuth  profile,  (c)  the  spectrum  estimated  using  enhanced 
SPIRE  and  (d)  azimuth  profile,  and  (e)  the  spectrum  estimated  using  MV  and  (f)  azimuth  profile. 
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5.  SPIRE  PERFORMANCE 


The  previous  section  dealt  with  computational  issues  relating  to  the  SPIRE  algorithm. 
In  this  section,  the  effect  of  various  uncontrolled  parameters  on  the  performance  of  the 
SPIRE  algorithm  is  studied.  These  parameters  are  discussed  in  the  next  few  sections 
and  include  SNR,  signal  spreading,  and  angular  spacing  between  signals.  Where 
appropriate,  the  results  using  the  MV  and  MUSIC  DF  algorithms  are  also  shown  for 
comparison  purposes. 


Since  the  SPIRE  algorithm  is  designed  to  estimate  the  direction  of  signal  power 
regardless  of  the  spatial  spreading  of  the  signal  source  (unlike  other  algorithms  such  as 
MUSIC  and  MV  which  assume  point  sources  only),  the  peaks  in  the  estimated  spatial 
spectrum  do  not  necessarily  yield  the  true  direction  of  spread-source  signals; 
particularly  when  the  ambiguity  problem  is  taken  into  account  (consider  Figure  1  If  for 
example).  A  more  appropriate  method  is  to  take  the  power  weighted  average  of  all  the 
grid  signals  associated  with  a  given  region.  For  example 

1  Mk 

(47)  (ftk  —  ^  v  Sik^ik 

i—  1 
1  Mk 

(48)  'Ipk  =  ^  v  Sik'fiik 

2=1 

where  <fik  and  ^  are  the  azimuth  and  elevation  bearing  estimates  of  the  kth  signal 
region,  Mk  is  the  number  of  grid  signals  required  to  model  the  region,  Sik  is  the  signal 
power  of  the  ith  grid  signal  in  the  kth  signal  region,  and  fak  and  ipik  are  the  associated 
azimuth  and  elevation  bearings  of  the  same  grid  signal. 

Although  no  research  has  been  done  to  determine  whether  the  above  bearing  estimation 
method  is  optimum,  for  simulation  testing  it  was  found  to  yield  results  which  were 
better  than  those  obtained  when  using  the  spectral  peaks. 


Signal  power  was  estimated  as  the  sum  of  the  grid  signals  in  the  associated  region,  that 
is 

Mk 

(49)  pk  =  ^  ^  S(k 

2=1 


where  pk  is  the  signal  power. 


In  this  section,  and  throughout  the  rest  of  this  report,  the  grid  size  used  for  the  SPIRE 
model  was  1°  x  1°,  unless  otherwise  specified. 

For  assessment  purposes,  the  processed  results  were  quantified  in  two  ways.  The  first 
was  the  measurement  of  the  failure  rate  of  signal  bearing  estimates,  and  the  second  was 
the  measurement  of  the  accuracy  of  the  successful  estimates.  A  bearing  estimate  was 
considered  to  be  a  failure  if  it  deviated  from  the  true  signal  bearing  by  more  than  half 
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the  array  beamwidth  (taking  into  account  both  the  azimuth  and  elevation  beamwidths). 
Accuracy  was  determined  by  calculating  the  root-mean-squared  (RMS)  error  of  the 
estimates  according  to 


(50)  RMS  Error  — 

where  the  summation  was  performed  for  all  H  successful  estimates  of  signal  m. 


Jj  -  4>m)2  +  {4>m  ~  VVn)2 


5.1  Effect  of  Noise 

The  effect  of  noise  was  investigated  through  simulation.  In  the  first  series  of 
simulations,  a  single  point-source  signal  at  (q p ,  ip)  =  (180°,  30°)  was  generated  and  the 
SNR  was  varied  from  -20  dB  to  +40  dB  in  2  dB  increments.  One  hundred  trials  were 
repeated  for  each  SNR  setting.  In  the  second  series  of  simulations,  the  point-source 
signal  was  changed  to  a  spread-source  signal  with  spread  parameters 
(A^,  Aip)  =  (30°,  15°),  but  all  other  parameters  remained  the  same.  The  results  after 
processing  the  data  with  the  SPIRE  algorithm  are  shown  in  Figure  16.  For  comparative 
purposes  the  results  using  the  SML  algorithm  are  also  shown. 

Several  features  of  the  results  are  worth  pointing  out.  The  SNR  at  which  the  failure  rate 
dramatically  departs  from  0%  is  called  the  threshold  SNR.  For  the  point-source  this 
occurs  between  -12  and  -10  dB  and  for  the  spread-source,  it  is  worse,  occurring  about 
-6  dB  for  the  SPIRE  algorithm  and  -8  dB  for  the  SML  algorithm.  The  poorer 
performance  for  the  spread-source  signal  (including  both  threshold  and  accuracy) 
compared  to  the  point-source  signal  is  a  function  of  the  amount  of  spreading. 

The  effect  of  signal  spreading  on  accuracy  can  be  attributed  to  two  factors:  signal 
model  uncertainty  and  the  filter  effect  (see  also  [1]).  The  first  factor,  signal  model 
uncertainty,  arises  from  the  fact  that  the  spatial  model  for  the  spread-source  signals  is 
stochastic  and  requires  a  sufficient  number  of  snapshots  to  build  up  the  appropriate 
statistics  in  the  data  as  was  discussed  in  Section  4.4.  For  a  single  point-source  signal,  a 
single  snapshot  is  sufficient.  For  a  spread-source,  the  number  of  snapshots  required  to 
achieve  a  given  accuracy  rises  as  the  spread  region  increases  in  size.  Conversely,  for  a 
given  number  of  snapshots,  accuracy  degrades  as  the  spread  region  increases. 

The  second  factor,  the  filter  effect,  can  be  understood  by  considering  that  many 
advanced  DF  algorithms,  such  as  MV,  MUSIC,  SML,  and  SPIRE,  can  be  interpreted  as 
techniques  which  work  by  designing  spatial  filters  to  reject  the  signal  content  of  the 
data.  Point-source  signals  are  matched  by  very  narrow  notch  filters  while  spread-source 
signals  are  matched  by  appropriately  shaped  band  rejection  filters.  The  greater  the 
amount  of  noise  rejected  by  the  filter,  the  greater  the  effect  on  the  estimation  error  since 
the  idea  is  to  reject  the  signal  but  pass  the  noise.  Hence  estimation  accuracy  degrades  as 
the  spread  region  increases. 

The  scales  used  to  display  the  RMS  bearing  errors  in  Figure  16(b)  were  chosen  because 
they  linearize  the  accuracy  results  for  the  point-source  signal  above  the  threshold  SNR. 
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gure  16:  Effect  of  the 
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In  this  region,  for  every  20  dB  increase  in  SNR  the  RMS  bearing  error  is  reduced  by  a 
factor  of  10.  Written  mathematically,  the  relationship  is  expressed  as 


(51)  eoc  VSNR 

where  e  is  the  RMS  bearing  error.  The  spread-source  also  exhibits  the  same  behaviour 
between  -8  and  0  dB,  but  above  this  SNR  the  error  begins  to  level  out  as  the  uncertainty 
in  the  signal  model  begins  to  dominate  the  error  (as  discussed  previously  in  Section 
4.4).  At  higher  SNR’s,  the  only  way  to  improve  accuracy  would  be  to  use  a  larger 
sample  size. 

Comparatively  speaking,  the  performance  of  the  SPIRE  algorithm  is  worse  than  the 
SML  algorithm.  For  the  point-source  signal,  the  failure  rate  was  the  same  for  both 
signals  as  was  the  accuracy  for  an  SNR  less  than  10  dB.  For  an  SNR  between  10  dB 
and  30  dB,  the  RMS  bearing  error  was  about  1.4  times  greater  for  the  SPIRE  algorithm 
than  the  SML  algorithm.  Above  30  dB,  the  increase  in  the  RMS  error  for  the  SPIRE 
algorithm  is  due  to  the  signal  model  grid  spacing  chosen  (1°  x  1°).  For  a  smaller  grid 
spacing  (i.e.  0.5°  x  0.5°),  the  SPIRE  error  reduces  from  12  times  to  1.4  times  the  SML 
error. 

For  the  spread-source  signal,  the  SPIRE  algorithm  again  produces  less  accurate  results 
with  RMS  errors  up  to  3  times  that  of  the  SML  algorithm. 

Although  the  poorer  performance  of  the  SPIRE  algorithm  might  at  first  seem 
disappointing,  it  is  a  consequence  of  the  greater  adaptability  of  the  SPIRE  signal  model. 
The  SML  algorithm  assumes  a  fixed  shape  and  power  density  for  each  signal’s  spatial 
profile.  In  these  examples,  this  profile  was  identical  to  what  was  actually  simulated 
giving  the  SML  algorithm  an  advantage.  In  the  real  world,  this  would  not  be  true  in 
many  situations  giving  the  SPIRE  algorithm  an  advantage  since  it  is  able  to  adapt  to  the 
conditions. 

Note  that  when  used  to  process  the  same  data,  the  MUSIC  and  MV  algorithms 
produced  the  same  results  as  the  SML  algorithm  for  the  point-source  signal,  but  were 
about  20 x  worse  for  the  spread-source  signal. 

The  last  set  of  results  displayed  for  the  noise  simulations  is  the  effect  of  signal-to-noise 
ratio  on  the  model  error  which  is  shown  in  Figure  17.  Comparing  the  two  signals, 
model  error  for  the  spread-source  is  relatively  independent  of  the  noise  level,  while  the 
model  error  for  the  point-source  increases  with  decreasing  noise  level.  This  is  due  to 
the  fact  that  the  point-source  cannot  be  perfectly  modeled  unless  a  signal  in  the  model 
grid  is  exactly  aligned  in  bearing  with  the  actual  signal  (in  the  simulations,  the  nearest 
grid  signal  was  misaligned  from  the  actual  signal  by  0.5°  in  both  azimuth  and 
elevation).  As  the  SNR  is  increased,  the  contribution  of  the  noise  component  decreases 
making  the  modeling  mismatch  in  the  signal  component  more  pronounced.  This  error 
could  be  reduced  by  decreasing  the  grid  size. 
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Figure  1 7:  The  effect  of  noise  on  model  error.  The  solid  line  shows  the  results  for  the  spread-source  signal  and  the 
dashed  lines  shows  the  results  for  the  point-source  signal. 


5.2  Detection  of  a  Weaker  Signal 

One  important  test  of  a  DF  estimator  designed  for  high  latitude  HF  operation  is  the 
ability  to  detect  a  weaker  point-source  signal  in  the  presence  of  stronger  spread-source 
signals.  To  evaluate  the  performance  of  the  SPIRE  algorithm,  simulations  were  run 
involving  a  single  spread-source  with  a  fixed  bearing  and  a  single  point-source  whose 
bearing  was  adjusted  incrementally,  beginning  with  a  large  initial  angular  difference, 
until  the  two  bearings  coincided.  After  each  increment,  the  signal  power  of  the 
point-source  was  increased  from  a  low  value  in  0.5  dB  intervals  until  the  failure  rate 
dropped  below  5%  (5  out  of  100  trials).  The  corresponding  SNR  of  the  point-source 
signal  is  defined  here  as  the  threshold  SNR  and  provides  a  good  indication  of  the  limits 
of  detectability  of  the  point-source  signal  for  the  given  signal  environment.  The 
relevant  signal  and  noise  parameters  are  shown  in  Table  1. 


Table  1:  Signal  Parameters  for  Signal  Detectability  Simulation 
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One  practical  difficulty  with  the  SPIRE  algorithm  is  that  as  the  separation  between  the 
two  sources  becomes  smaller  and  smaller,  the  regions  begin  to  merge  so  that  one  region 
describes  two  sources.  For  the  sake  of  these  simulations,  when  the  regions  had  merged, 
the  peak  in  the  combined  region  most  closely  corresponding  to  the  weak  signal  was 
used  to  determine  the  bearing  estimate  of  the  weak  signal.  This  was  considered 
reasonable  for  separations  of  15°  to  10°,  since  the  weaker  signal  was  clearly 
identifiable  when  examining  the  spatial  spectrum.  At  small  angular  separations, 
however,  no  attempt  was  made  since  separating  the  peak  of  the  weaker  signal  from  the 
peaks  associated  with  the  stronger  signal  was  essentially  a  hopeless  task. 
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The  results  from  the  simulations  are  shown  in  Figure  1 8.  In  this  case,  the  failure  rate  is 
not  shown  since  it  was  fixed  to  approximately  5%.  From  the  results  it  is  clear  that  when 
the  point-source  and  spread-source  signals  are  well  separated,  25°  or  more 
(>  1.2 4>bw),  the  presence  of  the  spread-source  signal  has  a  small  but  significant  effect. 
For  example,  the  threshold  SNR  for  a  point-source  without  any  other  signals  present 
occurs  at  -12  dB  (see  Figure  16).  In  the  presence  of  the  spread-source,  the  threshold 
SNR  ranges  from  -10  to  -8  dB,  a  degradation  of  2  to  4  dB.  As  the  separation  between 
signals  is  reduced  from  25°  to  10°,  the  threshold  SNR  increases  dramatically  to  8  dB 
which  is  20  dB  higher  than  when  no  spread-source  is  present. 

The  differences  between  the  wide  and  narrow  separation  cases  can  be  understood  in 
terms  of  using  a  spatial  filter  to  suppress  the  effects  of  the  spread-source  signal.  In  the 
wide  separation  case,  the  filtering  can  be  accomplished  easily,  leaving  only  the  noise  as 
the  main  source  of  error.  In  the  narrow  separation  case,  it  is  more  and  more  difficult  to 
filter  out  the  spread-source  signal  independently  of  the  point-source  signal  as  the 
separation  decreases.  Consequently,  the  spread-source  signal  begins  to  act  as  a  strong 
noise  background  and  the  threshold  increases  accordingly. 

The  accuracy  results  for  the  two  signals  are  relatively  constant  for  wide  spacing  and 
relatively  independent  of  the  presence  of  each  other.  For  example,  in  the  single  signal 
case,  the  accuracies  were  measured  to  be  0.33°  for  the  spread-source  signal  (20  dB 
SNR)  and  no  point  source  signal,  and  0.7°  for  the  point-source  signal  (-12  dB  SNR)  and 
no  spread-source  signal  -  worse  for  the  point-source  signal  due  to  the  much  lower  SNR. 
For  narrower  spacings,  the  threshold  signal  power  of  the  point-source  signal  increases, 
degrading  the  accuracy  of  the  spread-source  signal  until  the  separation  becomes  as 
small  as  8°.  The  accuracy  of  the  point-source  remains  in  the  range  of  0.5°  -  1.5°  for 
most  separation  angles  with  the  largest  errors  occurring  for  the  smaller  separations. 

The  simulations  were  also  repeated  using  both  the  SML  and  MUSIC  algorithms.  The 
assumed  number  of  signal  directions  for  the  SML  algorithm  was  two.  For  the  MUSIC 
algorithms,  six  directions  were  assumed,  since  occasionally  up  to  five  directions  were 
required  to  describe  the  spread-source  leaving  at  least  one  direction  for  the 
point-source.  In  most  cases,  however,  only  three  or  four  signal  directions  were  required 
for  the  spread-source  signal,  resulting  in  one  or  more  false  direction  estimates. 

Generally,  the  necessity  of  using  several  signal  directions  to  describe  a  spread-source 
signal  and  the  problem  of  extraneous  signals,  makes  the  interpretation  of  the  MUSIC 
results  somewhat  problematic.  For  the  sake  of  this  report,  only  the  two  (out  of  six) 
estimated  signal  bearings  closest  to  the  true  signal  bearings  were  used  when  generating 
the  statistical  results.  Additionally,  since  for  separation  angles  less  than  10°  there  was 
no  obvious  way  to  determine  whether  a  peak  was  associated  with  the  spread-source 
signal  or  the  point-source  signal  (like  the  SPIRE  algorithm),  no  statistical  results  were 
calculated  for  these  angular  separations. 

The  failure  and  accuracy  results  for  SML  are  shown  by  the  dashed  lines  in  Figure  18, 
and  the  corresponding  results  for  MUSIC  are  shown  by  the  dash-dot  lines. 
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(a) 


(b) 


(c) 

Figure  18:  Effect  of  angular  spacing  on  the  ability  to  detect  a  weaker  point-source  signal  in  the  presence  of  a 
stronger  spread-source  signal  showing  (a)  the  detection  threshold  SNR  for  the  weaker  point-source  signal,  (b) 
accuracy  of  the  spread-source  estimates  at  threshold,  and  (c)  accuracy  of  the  point-source  estimates  at  threshold. 
The  solid  lines  represent  the  SPIRE  results,  the  dashed  lines  represent  the  SML  results, and  the  dash-dot  lines 
represent  the  MUSIC  results. 
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(a) 


(b) 

Figure  19:  Detection  of  a  weaker  point-source  signal  in  the  presence  of  a  stronger  spread-source  signal  using 
MUSIC.  The  azimuth  spectrum  is  shown  for  an  elevation  angle  of  ip  =  30°.  The  simulation  parameters  for  the  noise 
and  spread-source  are  listed  in  Table  1,  and  the  point-source  parameters  were  (a)  <p  —  150°  with  a  signal  power  of 
-7.5  dB  and  (b)  <f>  —  172°  with  a  signal  power  of  14.5  dB.  When  the  separation  between  the  spread-source  and 
point-source  is  too  narrow,  as  in  (b),  there  is  confusion  as  to  which  spectral  peak  belongs  to  which  signal. 
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The  results  show  that  over  the  range  of  separation  angles  tested,  SPIRE  is  able  to  detect 
signals  2.5  dB  weaker  than  MUSIC.  Compared  to  SML,  performance  is  the  same  for 
separations  greater  than  20°,  but  generally  worse  for  smaller  separations.  Below  10°, 
the  SML  algorithm  is  the  only  algorithm  capable  of  separating  the  two  sources. 

The  accuracy  in  estimating  the  point-source  signal  direction  was  the  same  for  all  three 
algorithms.  However,  since  the  SML  and  SPIRE  algorithms  had  lower  thresholds,  this 
implies  that  tested  at  the  same  SNR  (e.g.,  the  MUSIC  threshold)  the  SML  and  SPIRE 
algorithms  are  more  accurate  than  MUSIC.  The  accuracy  of  MUSIC  for  estimating  the 
direction  of  the  spread-source  signal  was  very  poor,  highlighting  the  difficulty  of 
estimating  spread-source  signals  using  a  point-source  model.  The  accuracy  of  the 
SPIRE  algorithm  was  poorer  than  the  SML  algorithm,  indicating  that  some  degradation 
occurs  when  the  shape  of  the  signal’s  spatial  power  spectrum  is  not  known  a  priori. 

Generally  the  results  show  that  for  signal  environments  with  spread-source  signals, 
improved  modeling  leads  to  significantly  better  performance. 
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6.  PERFORMANCE  USING  HIGH  LATITUDE  OFF-AIR 
DATA  _ 


The  ultimate  test  for  any  algorithm  is  against  the  actual  data  for  which  the  algorithm 
was  designed.  In  the  following  sections,  an  analysis  of  the  SPIRE  algorithm  is  carried 
out  using  high  latitude  off-air  HF  data.  This  analysis  also  includes  results  from  the  MV 
and  MUSIC  algorithms. 

In  the  comparative  analysis  of  these  algorithms,  the  SPIRE  results  appear  to  be  worse 
than  the  MV  algorithm.  However,  it  was  found  that  mutual  coupling  problems  almost 
certainly  affected  the  data  measurements  possibly  making  the  results  for  the  MV 
algorithm  look  better  than  they  should  have  been.  The  diagnosis  of  this  problem  using 
the  SPIRE  error  function,  and  subsequent  discussion,  is  also  included  in  the  following 
analysis. 


6.1  The  Equipment 

The  measurement  system  used  to  collect  the  data  was  an  experimental  12  channel 
receiver  system,  called  “Vortex”,  located  at  CFB  Alert  on  the  northern  tip  of  Ellesmere 
Island  in  Northern  Canada  (82.50°  N,  62.35°  W).  The  receiver  system  was  connected 
to  an  antenna  array  which  had  the  geometric  configuration  shown  in  Figure  20.  The 
array  utilized  8  elevated  feed  monopole  antennas  from  the  inner  ring  of  a  Pusher  array 
(a  circular  array  with  24  antennas)  and  4  outlying  antennas  (also  elevated  feed 
monopole  antennas)  which  were  added  to  increase  the  array  aperture. 


Figure  20:  Three  dimensional  view  of  the  Vortex  antenna  array.  Each  grid  square  has  a  dimension  of  1 A  x  1  A. 
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The  Vortex  receivers  were  used  to  downconvert  the  input  signals  from  HF  to  2.5  kHz 
with  a  filtered  bandwidth  of  3.5  kHz.  The  downconverted  signals  were  then  digitized  at 
a  rate  of  10  kHz  and  the  data  stored  for  later  processing. 


6.2  Processing  Considerations 

To  generate  the  covariance  estimates,  an  FFT  was  performed  on  each  sample  block  of 
12  x  16000  data  points  and  only  the  positive  frequency  data  from  3.8  to  4.2  kHz  was 
retained.  This  served  the  dual  purpose  of  converting  the  data  to  IQ  format  and 
suppressing  interference  due  to  noise  and  other  unintended  HF  signals.  The  data 
covariance  matrix  was  then  formed  directly  from  10  consecutive  samples  of  this 
frequency  domain  data. 

Although  a  value  of  K  =  10  x  16000  (representing  1 6  seconds)  would  appear  to  be 
sufficient  based  on  the  discussion  in  Section  4.4,  this  does  not  take  into  consideration 
the  requirement  for  sample  to  sample  decorrelation.  The  time  required  to  achieve  this 
decorrelation  is  related  to  the  Doppler  spreading  of  the  signal,  and  can  be  approximated 
by  r  =  1  /{Doppler  Spread).  Using  a  Doppler  spread  value  of  40  Hz  (see  discussion  in 
Section  2.3),  this  gives  a  value  of  r  =  25  ms.  For  16  seconds  of  data,  this  corresponds 
to  an  effective  sample  size  of  K  =  640.  Since  40  Hz  represented  one  of  the  larger 
values  of  Doppler  spreading  observed  for  the  data  discussed  in  [10],  the  effective 
sample  size  for  the  data  discussed  here  may  have  been  lower.  For  this  reason,  the  1 6 
second  choice  for  the  sample  size  is  analyzed  in  more  detail  in  Section  6.5. 1 . 

When  processing  the  data  using  the  SPIRE  algorithm,  it  was  found  that  the  bearings  of 
peaks  in  the  estimated  spatial  spectrum  were  often  noticeably  different  than  bearings 
computed  using  the  weighted  region  approach  as  defined  in  (47)  and  (48).  This  was  due 
to  the  fact  that  the  estimated  signal  regions  were  often  asymmetric  and  sometimes 
contained  multiple  peaks.  In  some  cases,  it  also  appeared  that  the  peak  derived  bearings 
were  more  accurate  than  the  weighted  region  bearings.  Consequently,  peak  derived 
bearings  were  also  included  in  the  analysis. 

To  overcome  the  problem  of  quantization  for  the  peak  derived  bearings  (since  peak 
locations  were  restricted  to  bearings  of  the  signal  model  grid),  more  accurate  bearings 
were  interpolated  by  fitting  a  Gaussian  shape  to  the  grid  peak  and  its  immediate 
neighbouring  points.  The  peak  of  the  fitted  Gaussian  function  was  then  used  as  the 
improved  bearing  estimate. 

To  denote  the  difference  between  the  SPIRE  bearing  estimates  determined  using  the 
weighted  region  approach  and  bearing  estimates  determined  using  the  peaks  approach, 
the  former  are  termed  “region  bearings”  while  the  latter  are  termed  “signal  bearings” 
throughout  the  rest  of  this  report. 

For  processing  using  the  MUSIC  algorithm,  the  'number  of  signals’  parameter  was 
calculated  using  the  Akaike  Information  Criterion  as  described  in  [17]. 
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Both  the  MUSIC  and  MV  algorithms  were  adapted  to  use  the  measured  noise 
covariance  Cv  by  modifying  the  steering  vector  and  data  covariance  according  to 

C~*e 

(52)  e->  - fr - 

I!  <Ve  || 

(53)  R->  C^RC^ 

where  Cri  was  the  measured  noise  covariance  matrix. 

The  noise  covariance  Cr/  was  determined  in  the  same  way  as  the  signal  covariance 
matrix  except  all  the  positive  frequency  data  in  the  passband  was  retained  except  the 
signal  portion  from  3.8  to  4  kHz.  The  measurements  over  a  period  of  approximately  20 
minutes  were  then  averaged  to  ensure  a  stable  estimate 

6.3  The  Data 

Two  data  sets  were  chosen  for  this  discussion.  Both  collections  were  of  signals 
originating  from  the  CFH  transmitter  located  in  Halifax,  Nova  Scotia,  Canada.  The 
great  circle  signal  bearing  of  the  transmitter  was  181.9°.  The  received  signal  was 
collected  at  a  frequency  of  1 0.9445  MHz. 

The  first  data  set  was  collected  on  September  2,  1995  from  13:25:28  to  17:40:12  UT 
during  the  Arctic  daytime.  The  second  data  set  was  collected  on  January  24,  1996  from 
21 :01 :92  to  23:32: 12  during  the  Arctic  night  (which  lasts  24  hours  a  day  at  that  time  of 
year). 

Of  the  two  signal  periods,  the  first  represents  a  time  when  benign  signal  propagation 
conditions  would  normally  be  expected  (daytime)  while  the  second  represents  a  time 
when  disturbed  conditions  might  be  expected  [18].  This  is,  in  fact,  what  was  observed. 

6.4  DF  Results 

The  DF  results  for  the  two  data  sets,  with  the  three  different  algorithms,  are  shown  in 
Figures  21  and  22.  Additional  results  are  also  shown  for  the  SPIRE  algorithm  in 
Figures  23  and  24,  including  region  bearing  estimates  and  power,  for  both  data  sets.  For 
the  SPIRE  region  estimates,  results  were  discarded  if  the  corresponding  region  power 
was  less  than  0  dB. 

For  the  September  1995  data,  the  results  in  Figure  21  and  Figure  23  show  that  the  MV 
algorithm  produced  the  most  accurate  results  and  the  SPIRE  algorithm  produced  the 
worst  -  the  exact  opposite  from  what  was  expected. 

The  two  signal  bearings  estimated  by  the  SPIRE  algorithm  in  Figure  21a  and  b,  are 
peaks  from  the  same  region  so  that  the  scatter  in  the  estimates  is  a  reflection  of  the 
width  of  the  corresponding  region,  i.e.,  tens  of  degrees  of  in  both  azimuth  and  elevation. 
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The  large  amount  of  signal  spreading  explains  the  poorer  performance  of  the  MUSIC 
algorithm  compared  to  the  MV  algorithm.  MUSIC’S  superior  resolving  power  for 
point-source  signals  works  against  it  when  dealing  with  spread-source  signals.  The 
poorer  accuracy  of  the  SPIRE  algorithm,  however,  appears  to  be  due  to  entirely 
different  reasons  which  are  discussed  in  the  next  section. 

For  the  January  1996  data,  the  results  are  poor  for  all  three  algorithms  although  MV 
and  SPIRE  appear  to  produce  “cleaner”  results  in  Figures  22  and  24,  i.e.  appear  to  track 
moving  features  with  less  scatter.  The  greater  scatter  in  the  MUSIC  estimates  is 
consistent  with  idea  that  MUSIC  is  too  finely  tuned  for  point-source  signals. 

Comparing  the  bearing  estimates  to  the  power  levels  shown  in  Figure  24c,  the  times  of 
greatest  bearing  errors  occurred  when  the  signal  power  was  low  (less  than  0  dB) 
indicating  the  utility  of  making  these  measurements  as  a  means  of  qualifying  the 
measurements. 

Despite  the  worse  than  expected  bearing  accuracy,  the  SPIRE  results  do  provide 
valuable  insights  into  the  signal  environment.  For  example,  from  the  power  estimates 
for  September  1995  data  set,  there  appeared  to  be  only  one  dominant  signal  scattering 
region  active;  generally  in  the  great  circle  direction.  This  is  in  keeping  with  the 
expectation  that  the  signal  propagation  condition  would  be  benign  for  that  given  time 
period.  However,  comparing  the  region  bearing  estimates  in  Figure  23a  and  b  to  the 
signal  bearing  estimates  in  Figure  21a  and  b,  indicates  that  the  dominant  region  had  a 
complex  spatial  power  spectral  shape  (i.e.  more  than  one  peak)  and  that  it  was  quite 
dynamic  (i.e.  the  peaks  changed  location  from  estimate  to  estimate,  especially  the 
weaker  peak). 

Although  the  region  bearing  estimates  in  Figures  23  and  24  exhibit  less  scatter  than  the 
corresponding  signal  bearing  estimates  in  Figures  21  and  22,  the  accuracy  of  the  region 
bearing  estimates  appears  to  be  worse  overall.  The  weighted  approach  to  estimating 
region  bearings  was  based  on  the  premise  of  a  simple  hill-like  shape  for  the  spatial 
power  distribution  of  the  scattering  region.  The  complex  shaped  regions  and  dynamic 
nature  of  these  regions  may  make  the  weighted  approach  too  simplistic.  Better 
accuracy  may  entail  breaking  complex  shaped  regions  into  subregions  and  then 
estimating  the  weighted  bearing  of  each  of  these  subregions. 

Alternatively,  the  apparently  poorer  accuracy  of  the  weighted  region  approach 
compared  to  the  peak  approach,  and  the  apparently  poorer  accuracy  of  the  SPIRE 
algorithm  compared  to  the  MV  algorithm,  may  be  a  consequence  of  the  effect  of 
mutual  coupling  as  discussed  in  the  next  section. 

6.5  Problems  with  the  Data 

The  apparently  poorer  than  expected  performance  of  the  SPIRE  algorithm,  compared 
with  the  MV  algorithm,  suggests  a  possible  problem  with  either  the  SPIRE  signal 
model  or  the  way  the  data  was  generated  (equipment  problems).  To  address  this  issue. 
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Figure  21:  Azimuth  and  elevation  bearing  results  for  HF  high  latitude  off-air  data  collected  on  September  2 ,  1995 
using  the  SPIRE ,  MV  algorithm  and  MUSIC  algorithms :  The  azimuth  bearing  estimates  are  shown  in  (a)  for  SPIRE, 
(c)  for  MV,  and  (e)  for  MUSIC.  The  corresponding  elevation  bearing  estimates  are  shown  in  (b),  (d),  and  (f), 
respectively.  Only  the  two  strongest  signals  are  shown  with  black  representing  the  stronger  of  the  two  signals  and 
red  representing  the  weaker  of  the  two.  The  solid  horizontal  lines  in  (a),  (c),  and  (e)  represent  the  great  circle 
bearing. 
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Figure  22:  Azimuth  and  elevation  bearing  results  for  HFhigh  latitude  off-air  data  collected  on  September  2,  1995 
using  the  SPIRE,  MV  algorithm  and  MUSIC  algorithms.  The  azimuth  bearing  estimates  are  shown  in  (a)  for  SPIRE, 
(c)  for  MV,  and  (e)  for  MUSIC.  The  corresponding  elevation  bearing  estimates  are  shown  in  (b),  (d),  and  (f), 
respectively.  Only  the  two  strongest  signals  are  shown  where  black  represents  the  stronger  of  the  two  signals  and 
red  represents  the  weaker  of  the  two.  The  solid  horizontal  lines  in  (a),  (c),  and  (e)  represent  the  great  circle  bearing. 
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Figure  23:  More  results  of  processing  the  January  24,  1996  data  using  the  SPIRE  algorithm  and  showing  the 
estimated  (a)  region  azimuth  bearings,  (b)  region  elevation  bearings,  and  (c)  region  signal-to-noise  levels.  Only  the 
four  strongest  regions  are  shown  and  they  are  colour  coded  black,  red,  orange,  and  yellow  with  black  representing 
the  strongest  region  and  yellow  the  weakest  region. 
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Figure  24:  More  results  of  processing  the  January  24,  1996  data  using  the  SPIRE  algorithm  and  showing  the 
estimated  (a)  regions  azimuth  bearings,  (b)  region  elevation  bearings,  and  (c)  region  signal-to-noise  levels.  Only  the 
four  strongest  regions  are  shown  and  they  are  colour  coded  black,  red,  orange,  and  yeiiow  with  black  representing 
the  strongest  region  and  yellow  the  weakest  region. 
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several  different  aspects  of  model  and  data  generation  were  investigated,  namely,  signal 
correlation  and  sample  size,  spacing  and  position  of  the  signal  model  grid,  and  mutual 
coupling. 

6.5.1  Testing  Signal  Correlation 

The  amount  of  data  used  to  generate  each  data  covariance  matrix  was  assumed 
to  be  sufficient  for  decorrelation  of  the  signal  across  each  region.  To  test  the 
validity  of  this  assumption,  varying  amounts  of  the  September  1995  data  were 
used  to  generate  the  data  covariance  estimates  and  then  the  SPIRE  error  level 
determined.  Repeating  this  100  times  for  each  sample  size,  the  averaged 
model  error  results  are  shown  in  Figure  25. 


Figure  25:  Model  error  as  a  function  of  sample  size.  The  upper  curve  (solid  line)  shows  the  results  when  the 
September  2,  1995  data  is  processed.  The  lower  curve  (dashed  line)  shows  the  simulated  results. 

To  provide  a  reference  for  comparison  purposes,  simulated  data  was  also 
generated  and  the  averaged  model  error  computed  in  the  same  way.  The 
simulated  data  was  generated  according  to  the  expression 

(54)  X  =  WSN!  +  WnN2 
where  the  spatial  filter  matrices  were  defined  as 

(55)  Ws  =  (C  —  ct2R„)  2 
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I 

W—  R  2 

n  —  AVn 


(56) 

and  Ni  and  N2  were  spatially  white  noise  matrices.  The  quantity  C  —  a2Rn) 
was  the  noise  free  signal  model  estimated  by  the  SPIRE  algorithm  from  an 
arbitrary  sample  of  the  actual  data  (sample  length  of  16  seconds).  The  time 
domain  data  in  the  rows  of  Ni  was  filtered  to  have  the  same  spectral 
characteristics  shown  in  Figure  26,  which  is  an  example  of  the  spectral 
characteristics  of  the  CFH  signal  when  a  single  tone  was  being  broadcast.  The 
time  domain  data  in  the  rows  of  N2  was  not  filtered  (white  noise  only). 


Frequency  (Hz) 


(a) 


(b) 


Figure  26:  Spectral  characteristics  of  the  CFH  signal  are  shown  in  (a)  for  an  example  time  period  when  a  single  tone 
was  being  transmitted.  The  spectral  characteristics  of  the  filter  used  to  replicate  the  signal  part  of  the  spectrum  are 
shown  in  (b).  The  X-axis  in  both  plots  represents  the  frequency  with  respect  to  the  center  of  the  receiver  passband. 

The  purpose  of  generating  the  simulated  data  in  the  above  way  was  to 
replicate,  as  closely  as  possible,  the  actual  signal  characteristics  without 
violating  the  assumptions  made  about  signal  correlation.  Hence  Ni  would 
have  a  much  longer  correlation  constant  than  N2  due  to  the  temporal  filtering. 
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The  model  error  results  for  the  simulated  data  are  shown  in  Figure  25. 
Comparing  the  actual  and  simulated  results  confirms  that  there  is  a  problem 
with  the  model,  but  suggests  that  it  was  probably  not  due  to  insufficient 
correlation.  For  sampling  times  greater  than  about  1  second,  the  modeling 
error  for  the  off-air  data  begins  to  flatten  out  while  the  modeling  error  for  the 
simulated  data  continues  to  decrease.  This  might  suggest  that  an  extremely 
long  correlation  time  was  involved,  however,  even  combining  data  from 
different  days/months  did  not  significantly  reduce  the  model  error  showing 
that  correlation  was  not  the  problem. 

6.5.2  Testing  the  Model  Grid 

Another  possible  explanation  is  that  the  recorded  signals  in  the  data  had  a 
strong  point-source  component  so  that  misalignment  of  the  model  grid 
contributed  to  the  error,  as  discussed  in  Section  5.1.  However,  reducing  the 
model  grid  spacing  to  0.5  x  0.5  had  no  effect  on  the  model  error,  nor  did 
adjusting  the  grid  up  to  ±0.5°  in  azimuth  and  elevation  -  essentially  ruling  out 
misalignment  as  the  source  of  model  error. 


6.5.3  Testing  for  Coupling  Effects 

Having  ruled  out  problems  which  could  be  mitigated  by  appropriately 
processing  the  data,  the  most  likely  problem  was  antenna  coupling  with  the 
local  environment  (local  site  multipath)  as  well  as  coupling  with  other 
antennas  in  the  array  (mutual  coupling).  Coupling  problems  were  extremely 
likely  given  that  the  local  terrain  was  neither  perfectly  flat  nor  the  ground 
perfectly  conducting  (frozen  ground  has  low  conductivity),  and  that  the 
antennas  in  the  inner  Pusher  ring  are  spaced  closer  than  1/4  wavelengths  at 
10.9445  MHz,  virtually  guaranteeing  antenna  to  antenna  coupling. 

To  investigate  further,  the  effects  of  antenna  coupling  were  added  to  the 
simulation  models  using  the  method  outlined  in  [19]  where  impedances  were 
calculated  using  the  approach  in  [20],  At  best,  this  can  only  be  considered  an 
approximation  since  the  antenna  is  modeled  to  be  an  ideal  monopole  on  a 
perfectly  conducting  ground.  In  reality,  the  antennas  were  elevated  feed 
monopoles  on  poorly  conducting  ground.  Hence,  the  main  purpose  of  the 
simulation  was  to  highlight  the  kinds  of  errors  that  multipath  causes,  not 
replicate  the  exact  errors  that  actually  occurred. 

Using  modeling  errors  as  a  guide,  the  initial  simulations  showed  that  the 
modeled  mutual  coupling  voltage  levels  were  too  high,  so  these  levels  were 
reduced  to  30%  of  their  unmodified  values.  Using  this  adjusted  mutual 
coupling  model,  the  theoretical  effect  of  sample  size  on  model  error  is  shown 
in  Figure  27.  In  this  case  there  is  much  better  agreement  between  the 
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simulated  results  and  the  actual  results  than  in  Figure  25  strongly  implicating 
antenna  coupling  as  the  source  of  errors. 


Figure  27:  Model  error  as  a  function  of  sample  size.  The  upper  curve  (solid  line)  shows  the  results  when  the 
September  2,  1995  data  is  processed.  The  lower  curve  (dashed  line)  shows  the  simulated  results  when  mutual 
coupling  effects  are  included. 

To  illustrate  how  modeling  errors  translate  into  bearing  errors,  an  example  of 
the  effects  of  mutual  coupling  on  the  SPIRE  estimated  power  spectrum  is 
shown  in  Figure  28b.  Comparing  this  to  the  ideal  spectrum  in  Figure  28a,  the 
major  effect  observed  is  the  distortion  of  the  shape  of  the  signal  region  and  the 
introduction  of  false  regions.  The  shape  distortion  leads  to  biasing  in  the 
bearing  estimates  while  the  false  regions  leads  to  false  estimates  (i.e.  wild 
bearings).  In  the  worse  case,  distortion  can  lead  to  extra  peaks  in  the  main 
signal  region  causing  potential  errors  for  peak  search  methods. 

Given  the  distorting  effects  of  antenna  coupling,  it  would  be  reasonable  to 
suggest  that  the  complex  nature  of  the  power  spectrum  shown  for  the  off-air 
example  Figure  28c  was  mainly  due  to  the  effect  of  multipath  -  the  true  power 
spectrum  being  simpler  in  nature.  Unfortunately,  the  evidence  is  not  strong 
enough  to  make  this  assessment  with  any  certainty. 

More  compelling  evidence  is  shown  in  Figure  29.  Five  hundred  samples  of 
data  were  generated  using  the  same  signal  model  used  in  Figure  28  (including 
the  mutual  coupling  effects)  except  that  the  azimuth  bearing  was  slowly 
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?  28:  The  effect  of  mutual  c 
\ted  SPIRE  spectrum  whet 
iple  from  the  September  1 


changed  from  171.9°  to  191.9°  in  a  uniform  manner.  Estimates  were  made 
using  the  SPIRE  algorithm  (weighted  and  peak  approach)  and  the  MV 
algorithm.  The  true  bearing  is  shown  as  the  solid  line  in  each  of  the  plots. 

The  most  remarkable  thing  about  these  results  is  that  for  samples  150  to  500 
the  MV  algorithm  produced  estimated  bearings  of  180°  ±  1°  even  though  the 
actual  bearings  changed  from  175.9°  to  191.9°.  This  “flattening”  effect  is  also 
apparent  for  the  SPIRE  algorithm  results  when  the  peak  method  was  used,  but 
not  as  severe.  When  the  weighted  region  method  was  used,  the  flattening 
effect  disappears.  There  is  some  biasing  of  up  to  0.5°,  but  this  is  substantially 
less  than  the  results  for  the  peak  method  or  the  MV  algorithm. 

The  main  implication  is  that  the  superior  performance  of  the  MV  algorithm 
for  the  off-air  data  (particularly  the  September  1995  data),  might  conceivably 
have  been  an  illusion.  That  is,  it  is  entirely  possible  that  the  great  circle 
bearing  was  coincidentally  favoured  due  to  the  effects  of  mutual  coupling 
making  the  MV  results  appear  better  than  they  should  have  been. 

Regardless  of  the  relative  performances  observed,  it  seems  highly  probable 
that  the  results  for  all  the  algorithms  were  adversely  affected  by  mutual 
coupling.  Hence,  a  more  accurate  assessment  of  algorithm  performance  is  not 
possible  unless  coupling  effects  can  be  calibrated.  Ideally,  this  could  have 
been  accomplished  using  an  airborne  transmitter  to  measure  the  antenna  array 
response  as  a  function  of  azimuth  and  elevation  angles  for  each  frequency  of 
interest.  Unfortunately,  due  to  logistical  and  cost  reasons,  this  wasn’t  done. 
Calibration  methods  based  on  the  data  itself  have  been  developed,  but  to  the 
author’s  knowledge,  none  are  appropriate  for  the  high  latitude  signal 
environment.  Hence  more  research  is  required  in  this  area. 
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(a) 


(b) 


(c) 

Figure  29:  The  effect  of  mutual  coupling  on  bearing  estimation  for  a  moving  spread-source  showing  the  azimuth 
bearings  estimates  for  (a)  the  SPIRE  algorithm  using  the  weighted  region  approach,  (b)  the  SPIRE  algorithm  using 
the  peak  approach,  and  (c)  the  MV  algorithm.  The  true  bearings  are  shown  by  the  solid  line. 
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7.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  Spatial  /ncoherent  Region  Estimator  (SPIRE)  is  a  new  algorithm  based  on 
maximum  likelihood  principles  which  was  developed  to  estimate  the  spatial  power 
spectrum  using  measurements  from  an  N-channel  antenna  array.  The  algorithm  is 
distinguished  from  other  spectral  estimation  algorithms  in  that  it  assumes  that  all 
signals  are  uncorrelated  but  makes  no  assumption  about  their  spatial  shape  and  power 
distribution. 

Several  user  controlled  aspects  of  the  algorithm  were  examined  through  simulation 
testing.  This  included  resolving  spatial  ambiguities,  determining  the  most  appropriate 
grid  spacing  for  the  signal  model  and  the  minimum  sampling  size,  and  exploring  the 
limit  on  the  maximum  number  of  signals  that  could  be  handled. 

The  testing  showed  that  the  easiest  way  to  deal  with  spatial  ambiguities  is  by  choosing 
the  simplest  signal  model  which  fits  the  data.  The  testing  also  showed  that  the  grid 
spacing  and  sampling  size  parameters  could  be  varied  for  a  range  of  values  over  which 
estimation  performance  was  either  virtually  unaffected  or  predictable.  Additionally,  it 
was  also  shown  that  the  SPIRE  algorithm  can  estimate  a  greater  number  of  signals  than 
there  are  antennas  in  the  antenna  array,  exceeding  the  traditional  N  —  1  limit  (where  N 
is  the  number  of  antennas).  However,  practically  speaking,  remaining  below  this  limit 
is  more  advisable. 

Simulation  testing  of  the  SPIRE  algorithm  was  also  carried  out  to  evaluate  its 
performance  as  a  function  of  various  environmental  conditions.  This  included  the  effect 
of  noise,  signal  spreading,  and  the  detection  of  a  weaker  signal  in  the  presence  of  a 
stronger  signal.  In  all  cases,  the  performance  of  the  SPIRE  algorithm  was  predictable. 
Accuracy  was  similar  to  other  superresolution  algorithms  when  dealing  with  signals 
with  no  spatial  spreading,  but  better  when  there  was  spreading. 

Finally,  testing  was  performed  using  off-air  data  collected  at  the  Arctic  site  CFS  Alert. 
The  results  were  inconclusive  as  antenna  mutual  coupling  effects  were  found  to  have 
corrupted  the  data.  The  analysis  of  the  data  did,  however,  showcase  the  advantages  of 
the  SPIRE  algorithm  in  helping  diagnose  the  non-ideal  nature  of  the  off-air  data. 

Generally  it  was  demonstrated  that  the  SPIRE  algorithm  is  able  to  estimate  the  spatial 
power  spectrum  of  the  radio  environment  with  a  higher  resolution  and  better  accuracy 
than  previously  possible.  As  an  analytical  tool,  the  SPIRE  algorithm  provides  a 
powerful  new  method  for  analyzing  the  spatial  nature  of  signals,  as  well  as  a  means  of 
detecting  data  and  equipment  problems.  Although  originally  developed  for  the  high 
latitude  HF  signal  environment,  the  algorithm  can  be  applied  to  any  TV-channel  data  set 
provided  the  signals  are  uncorrelated  both  temporally  and  spatially. 

Mutual  coupling  of  antennas  in  the  measurement  array  remains  a  problem.  Although 
there  is  reason  to  believe  that  the  SPIRE  algorithm  may  be  less  sensitive  to  mutual 
coupling  problems  than  other  algorithms,  it  is  still  adversely  affected.  For  the  high 
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latitude  HF  measurements  made  from  CFS  Alert  in  the  mid  1990’s,  calibrating  the 
coupling  effects  is  of  paramount  importance  if  the  potential  usefulness  of  this  data  set  is 
ever  to  be  realized.  Since  no  in  situ  calibration  was  ever  done,  calibration  would 
involve  deriving  the  correction  coefficients  from  the  data  itself.  Unfortunately  no  good 
method  has  been  developed  to  do  this  kind  of  calibration,  so  continued  research  in  this 
area  is  required.  Additionally,  since  in  situ  calibration  is  often  expensive  and  difficult, 
many  modem  DF  systems  would  also  benefit  greatly  from  this  kind  of  research. 
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